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FOREWORD 


System  433L;  project  1.0;  task  1.2.  This  TR  has  been  prepared  for  United 
Aircraft  Corporation,  East  Hartford,  Conn.,  under  Subcontract  No.  15107  to 
Contract  No.  AF19(628)-3437,  by  The  Travelers  Research  Center,  Inc., 

250  Constitution  Plaza.  Hartford,  Conn.  The  Research  Center’s  publication 
number  is  7463-182.  Robert  L.  Houghten,  Lt.  Colonel,  USAF,  is  Acting  System 
Program  Director.  This  report  covers  the  period  January— August,  1965,  and 
was  submitted  for  approval  on  October  6,  1965. 


ABSTRACT 


A  physical  model  for  the  prediction  of  large-scale  low  cloudiness  [ESD-TR-65-3] 
was  modified.  Results  of  three  12-hour  forecasts  (for  different  synoptic  situations), 
computed  with  the  new  model,  are  presented. 

Equations  for  the  contact  layer  were  derived  for  the  forced  convection,  free 
convection,  and  strong  inversion  regimes.  The  requirement  for  upper  boundary  condi¬ 
tions  is  reduced;  only  the  geostrophic  wind  components  and  upper-level  cloudiness  are 
required  from  a  free-air  model.  Empirical  methods  were  developed  for  computing 
the  instrument  shelter  level  temperature  and  relative  humidity.  The  eddy  diffusion 
coefficient  for  heat  and  water  is  a  function  of  the  Richardson  number  throughout  the 
boundary  layer.  The  linear  geostrophic  wind  shear  within  the  boundary  layer  is  com¬ 
puted  from  the  predicted  temperature. 

The  model  has  not  been  thoroughly  tested:  preliminary  results  suggest  that,  given 
a  careful  analysis  of  routine  observational  data,  this  boundary  layer  model  will  serve 
operationally-meaningful  diagnostic  and  predictive  functions. 

The  logic  used  in  constructing  the  computer  program  for  the  model  is 
presented. 


REVIEW  AND  APPROVAL 

Publication  of  this  technical  report  does  not  constitute  Air  Force  approval  of  the 
report’s  findings  or  conclusions.  It  is  published  only  for  the  exchange  and  stimulation 
of  ideas. 


Lt.  Colonel,  USAF 

Acting  System  Program  Director 
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SECTION  I 


INTRODUCTION 

It  is  frequently  noted  that  physical  atmospheric  prediction  models  rarely  have 
merited  designation  as  weather  prediction  models.  Empirical  investigations  have  shown 
that  the  interpretation  of  even  an  accurate  circulation  prediction  in  terms  of  weather 
conditions  is  not  a  trivial  problem. 

Some  ten  years  ago,  Smagorinsky  and  Collins  [26]  advanced  a  procedure  for 
utilizing  a  two-parameter,  quasi-geostrophic  model  for  predicting  precipitation  amounts. 
The  hypothesis  that  diabatic  effects  might  have  a  significant  role  in  cases  of  strong, 
baroclinic  development  led  many  investigators  to  develop  physical  prediction  models  that 
incorporated  the  condensation  process  more  or  less  along  the  lines  suggested  by 
Smagorinsky  and  Collins.  Gambo  [11]  has  indicated  recently  that  diabatic  effects,  while 
they  may  not  lead  to  rapid  deepening,  can  be  very  significant  in  the  maintenance  of  accurate 
phase  relations  in  multi-level  circulation  models.  At  the  present  time,  neither  the  Air 
Force  nor  the  National  Meteorological  Center  (NMC)  errplcys  physical  prediction  models 
that  incorporate  the  diabatic  processes  which  are  intimately  related  to  weather  phenomena. 
In  general  circulation  research,  the  utilization  of  diabatic  models  is  clearly  essential,  and 
in  the  light  of  the  requirement  for  long-range  weather  predictions,  operational  prediction 
will  have  to  be  reoriented  to  the  utilization  of  diabatic  physical  prediction  models. 

The  region  of  the  troposphere  which  is  most  greatly  influenced  by  diabatic  processes 
is  the  planetary  boundary  layer.  Because  most  horizontally-extensive  low  cloudiness  is 
confined  to  this  layer,  it  appears  desirable  to  approach  the  problem  of  the  prediction  of 
synoptic-scale  low  cloudiness  through  the  use  of  a  reasonably  complete  physical  model  of 
the  boundary  layer.  The  model  which  is  presented  in  this  paper  is  an  effort  in  this  direc¬ 
tion. 

Some  work  has  been  done  previously  in  the  formulation  of  subsynoptic-scale  physical 
models  of  the  boundary  layer,  notably  by  Fisher  [8]  and  Estoque  [7].  Their  work  could 
not  be  directly  extended  because  of  their  requirement  for  non-routine  meteorological 
observations.  It  was  necessary,  therefore,  to  consider  the  feasibility  of  constructing  a 
meaningful  model  which  would  utilize  only  routine  observational  data  but  yet  incorporate 
the  significant  physical  processes. 


1 


✓ 

A  number  of  suggestions  made  to  the  writer  by  Drs.  Arnason  and  Pandolfo  of  The 
Travelers  Research  Center,  Inc.  (TRC),  Prof.  Davidson  of  New  York  University,  and 
Prof.  Blackadar  of  Pennsylvania  State  University  were  helpful  in  indicating  the  extent  to 
which  certain  physical  approximations  might  be  justified.  Foremost  among  the  modeling 
simplifications  adopted  in  our  work  are: 


(a)  the  use  of  a  static -diagnostic  horizontal-wind  equation, 

(b)  the  use  of  an  empirical  surface  temperature  prediction  technique, 

(c)  the  use  of  empirical  formulas  for  the  stress  and  geostrophic 
deviation  angle, 

(d)  the  neglect  of  radiative  heat  flux  divergence, 

(e)  the  use  of  the  constant-flux  contact  layer, 

(f)  the  use  of  an  empirical  surface  humidity  prediction  technique. 


The  research  presented  in  this  paper  was  essentially  oriented  to  a  specific  applica¬ 
tion,  i.e.,  low-cloud  prediction.  A  broader  view  of  the  framework  within  which  this 
modeling  approach  may  be  applied  has  been  indicated  above,  but  has  not  yet  been  pursued. 

After  development  of  the  model  to  its  present  form,  an  extensive  period  of  careful 
testing  is  desirable.  From  a  practical  viewpoint,  these  tests  should  indicate  the  relative 
usefulness  of  the  model  in  comparison  with  other  methods  for  low-cloud  prediction.  An 
equally  useful  purpose  for  such  tests  would  be  to  unveil  the  hidden  characteristics  of  the 
model  atmosphere.  Within  the  scope  of  the  effort  reported  here,  we  have  been  able  to 
conduct  only  a  few  tests  of  the  model. 
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SECTION  U 


DERIVATION  OF  MODEL  EQUATIONS 


The  equations  governing  the  model  atmosphere  are  written  in  Cartesian  coordinates. 
The  horizontal  coordinates,  x  and  y,  are  rectangular  coordinates  on  a  polar  stereographic 
map  projection.  Because  the  pertinent  physical  boundary  is  the  surface  of  the  ground  or 
ocean,  and  not  a  fictitious  mean-sea-level,  we  chose  the  height  above  the  terrain  as  the 
vertical  coordinate,  Z. 

The  planetary  boundary  layer  is  decomposed  into  two  sublayers.  A  shallow  layer, 
denoted  as  the  contact  layer,  occupies  the  region  between  Z  =  0  and  Z  =  h  (50  m).  The 
bulk  of  the  boundary  layer  is  termed  the  transition  layer.  It  has  been  assumed  that  the 
upper  boundary  of  the  planetary  boundary  layer  may  be  set  at  2  km  above  the  terrain 
height. 

1.  Horizontal  Equations  of  Motion 

We  assume  that  the  horizontal  wind  components,  u  and  v,  are  governed  by  the 
balance  equations, 


(II-l) 


(II- 2) 


where  oi,  the  square  root  of  the  ratio  of  the  Coriolis  parameter  to  twice  the  eddy 

viscosity,  is  assumed  to  be  independent  of  Z,  and  u  and  v  ,  the  geostrophic  wind  com- 

g  g 

ponents,  are  assumed  to  be  linear  functions  of  Z.  Under  the  boundary  conditions, 


u  — -  u 


g 


as  Z  —  00 


(II- 3 ) 


g 


and 


u  —  U 


at  Z  =  h 


(II— 4) 


v  =  V 


the  solutions  of  Eqs.  (II-l)  and  (II- 2)  may  be  given  by 
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(U-5) 


u  =  u  +  e  ^  {(U  -  u^1)  cos  [a(Z-h)]  +  (V  -  v*1) 

s  s  s 

sin  [a(Z-h)]} 

-Q;(Z-h)  ,/ir  h  r  ...  h 

v  =  v  +  e  {(V  -  v  )  cos  [o;(Z-h)]  -  (U“  u  ) 

S  S  5 

sin  [a(Z-h)]} 


(II- 6) 


h  li 

in  which  u  and  v  are  the  values  of  u  and  v  at  Z  =  h.  The  values  of  U  and  V  required 
g  g  g  g 

in  Eqs.  (II— 5)  and  (II— 6)  are  obtained  from  formulas  derived  subsequently. 

2.  The  Geostrophic  Wind 


The  geostrophic  wind  is  assumed  to  be  a  linear  function  of  Z: 

u  =  uH  +  B  (H-Z) , 
g  g 

v  =  vH  +  C  (H-Z) . 
g  g 


(H-7) 


(H-8) 


Using  the  procedure  outlined  in  Haurwitz  [14],  one  may  relate  B  and  C  to  the  horizontal 
temperature  gradient  within  the  transition  layer.  Taking  appropriate  account  of  the 
terrain  variation,  we  may  derive  the  following  relations: 

H 


u  =  u  + 
g  g 


f  h  ay  IT 


(H-9) 


dZ}  , 


v  =  v 
g 


H  Ph  -  z“|  H 
■ +  ' 

,H  _8  |lL,, 

f  4  ax  1  t  J  ’ 


T  -  T 
h  H 


H 


(II- 10) 


in  which 


u  =  - 


qg  9E 
f  9y  ’ 


qg  9E 

V  =  — 

f  ax’ 


<n-ii) 

(11-12) 


H  H 

and  Ug,  v^  are  the  geostrophic  wind  components  at  Z  =  H,  T^  and  T^  are  the  tempera¬ 
tures  at  Z  =  h  and  Z  =  H,  E  is  the  height  of  the  terrain  above  mean  sea  level  and  q  is  the 
map  scale  factor. 
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3_. _ The  Continuity  Equation 

Neglecting  the  partial  derivative  of  air  density  in  the  continuity  equation  leads 
to  the  simplified  equation, 


9w 

dz 


=  -  <7 


du  3v 
9x  9y , 


(11-13) 


Since  w  must  vanish  at  Z  =  0  due  to  the  viscosity  of  the  air,  we  may,  to  good  approxima¬ 
tion,  solve  Eq.  (II-9)  using  the  boundary  condition,  w  =  0  at  Z  =  h. 


fZ  /  9u  3v  \  _  /TT  .  . . 

w“  ■  "M**  %}  (  ’ 

Due  to  our  choice  of  coordinate  system,  the  vertical  velocity,  w,  differs  from  the 
velocity  normal  to  surfaces  which  parallel  mean-sea-level.  If  the  latter  is  denoted  by  w, 
we  may  write, 

oj  -  w  +  w ,  (11-15) 


with 


(11-16) 


4.  The  Heat  Transfer  Equation 

The  first  law  of  thermodynamics,  the  ideal  gas  law,  and  the  hydrostatic  equation 
may  be  combined  to  derive  a  heat  transfer  equation, 


where 


dT  _  RT  d£ 

dt  C  p  dt 
P 


R 

P 

T 

K 


H 


g 

'  dT  J 
l  dt  ) 


i(KH 


S—  K 
C  T  H 


91  +  _g_ 
az  c 

p=! 

9T  +  _g_ 

az  c 


} 


IdT 

(  dt 


(11-17) 


w 


is  the  specific  heat  of  air  at  constant  pressure 

is  the  gas  constant  for  dry  air 

is  the  air  pressure 

is  the  air  temperature 

is  the  eddy  conductivity 

is  the  gravitational  constant 

is  the  heating  due  to  water  substance  phase  changes 


We  simplify  Eq.  (11-12)  by  the  approximation, 
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(w  +  w) 


(11-18) 


HT  d£  _  _  _g_ 

C  p  dt  C 
P  P 

and  by  neglecting  the  third  term  on  the  right  hand  side  of  Eq.  (11-17).  The  resulting 
equation  is 


(11-19) 


5.  The  Water  Substance  Transfer  Equation 

The  ratio  of  the  water  vapor  density  to  the  total  density  of  the  mixture  of  air,  water 
vapor,  and  liquid  water  is  defined  as  the  specific  humidity  and  denoted  by  q.  We  simil¬ 
arly  define  the  specific  moisture,  r,  as  the  ratio  of  water  substance  (vapor  and  liquid) 
density  to  the  density  of  the  mixture.  When  precipitation  is  omitted  from  consideration, 
we  may  write  the  following  equations  for  q  and  r. 


dr  _ 

f  3r 

in  -  i 

3r 

dt  ~  ° 

U  _  “T 

[_  3x 

v9y 

*--<r 

l>+ 

vfa 

at 

1  8* 

dy 

]- 

]■ 


ar 

w  —  + 

az 


9q  a 

W  — ^  +  - 

az  az 


_arK  -^1 

az  |_  v  az  J 

Ml-  isi 


(11-20) 


(11-21) 


where  is  the  eddy  diffusivity  for  water  vapor  and  liquid  water.  The  approximation 


that  both  vapor  and  droplets  follow  the  air  motion  is  convenient  and  in  accord  with 

l» 


observations  [21). 


represents  the  rate  of  conversion  of  vapor  to  liquid  in  con¬ 


version  of  vapor  to  liquid  in  condensation. 


6.  Condensation,  Evaporation,  and  Latent  Heat  Exchange 

In  the  heat  transfer  equation  and  the  water  vapor  transfer  equation,  source  terms 
appear  for  which  a  computational  procedure  is  required.  We  have  adopted  the  method 
used  by  Fisher  and  Caplan  [8]  modified  after  the  fashion  suggested  by  McDonald  [22]. 

In  our  technique,  we  use  a  standard  pressure— height  relationship  to  provide  the  air 
pressure  value  needed  in  the  computation  of  the  saturation  value  of  specific  humidity. 
This  relation  is  given  by 

—  c 

p  =  [1.013  -  1.065  X  10"  (Z  +  E)J  (11-22) 
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in  which  p  is  in  bars,  and  Z  and  E  are  in  cm.  The  saturation  specific  humidity, 
qs,  at  a  temperature,  T(°K),  and  at  a  height,  Z,  is  given  to  good  approximation  by 

3.8  X  10"3 

q  - - a - 

[1.013  -  1.065  X  10  (Z+E)] 


(11-23) 


7.  The  Coefficients  of  Eddy  Diffusivity  and  Eddy  Viscosity 

The  computational  formulas  used  to  compute  the  coefficients  of  eddy  diffusivity 
and  eddy  viscosity  within  the  transition  layer  are  extensions  of  formulas  which  have 
been  derived  from  similarity  theories  for  the  contact  layer.  Richardson’s  number, 
R.,  is  defined  as 


=  se  /  85,: 
e  sz  '  »»  1 


i  e  az.  '  0Z 

The  free  and  forced  regimes  of  turbulent  convection  are  discriminated  by  the  use  of 
a  critical  value  of  R^.  When  Rj  s  -  0.03  we  assume  that  free  convection  exists; 
when  Rj  >  -  0.03  we  assume  that  forced  convection  exists. 

The  eddy  diffusivities  for  heat  and  water  (both  vapor  and  liquid)  are  taken,  to 
be  equal.  If  free  convection  prevails,  we  write, 


(11-24) 


K  =  Kf ,  =  \Z 
v  H 


JL 

e  az 


1/2 


(11-25) 


If  forced  convection  prevails,  we  write, 


Ky  =  Kh  =  [kZ  (1-/3R ,)Y 


2  |  8? 
az 


(11-26) 


The  values  computed  from  these  formulas  are  adjusted  to  lie  between  10^  and  10^ 

cm2  sec  *.  If  R  >  0  and  I  ^  1  =  0,  or  (1-/3R.)  ^  0,  we  set  the  coefficients  to  the 
1  oZ  1 

minimum  value. 

The  coefficient  of  eddy  viscosity  is  assumed  to  be  invariant  with  height  and  equal 
to  the  value  computed  for  it  from  the  contact-layer  equations  at  Z  =  h. 


8.  Boundary  Conditions 

The  equations  presented  above  require  for  solution  the  specification  of  boundary 
and  initial  conditions.  The  initial  conditions  can  be  specified  from  observational  data 
throughout  the  region  of  integration.  Lateral  boundary  conditions  are  required  for  those 
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equations  involving  horizontal  advection  terms.  These  time-dependent  conditions  are 
needed  only  on  those  portions  of  the  boundary  through  which  the  horizontal  flow  is 
directed  into  the  region  of  integration.  Errors  in  the  specification  of  conditions  on 
these  boundaries  will  propagate  into  the  interior  of  the  region  at  approximately  the  speed 
of  the  air. 

Boundary  conditions  are  also  required  for  the  heat  and  water  transfer  equations  on 
the  vertical  boundaries.  On  the  lower  boundary  we  specify  the  eddy  flux  of  heat  and 
water  substance.  The  computation  of  these  boundary  values  is  carried  out  using  the 
formulas  derived  subsequently  through  consideration  of  the  properties  of  the  contact 
layer.  On  the  upper  boundary  q,  T  and  r  are  computed  from  simplified  forms  of  the 
basic  equations.  The  simplification  involves  the  neglect  at  this  boundary  of  the  convergence 
of  the  eddy  flux. 

The  vertical  boundary  conditions  needed  for  the  complete  solution  of  the  horizontal 

H  H 

wind  equations  are,  first,  the  specification  of  the  temporal  variation  of  u  and  v  ,  and 

S  S 

second,  the  speciffciioncf  U  and  V  at  Z=h.  The  first  conditions  may  be  obtained  from  fore¬ 
casts  of  the  geopotential  made  by  dynamical-prediction  models  at  appropriate  pressure 
levels.  U  and  V  are  computed  from  formulas  derived  using  the  contact  layer  equations. 

9.  Constant  Flux  Hypothesis 

The  rate  at  which  meteorological  properties  are  transferred  by  eddy  motions  within 
the  region  near  the  ground  is  so  large  (at  least  at  certain  times)  that  it  is  quite  clear 
from  standard  observations  that  this  eddy  transport  must  be  very  nearly  non-divergent. 

For  example,  the  eddy  heat  flux  has  been  measured  to  be  several  hundred  milli-langleys 
per  minute.  If  that  amount  of  heat  transfer  were  to  converge  within  a  layer  of  a  few  hundred 
feet  thickness,  it  would  produce  temperature  changes  throughout  the  layer  of  several  tens 
of  degrees  (centigrade)  in  an  hour.  Detailed  measurements  of  the  eddy  transport  of  heat, 
vapor,  and  momentum  have  been  made  utilizing  instruments  mounted  on  towers  at  various 
heights.  To  the  accuracy  obtainable  with  these  instruments,  it  has  been  found  (19]  that 
the  convergence  of  the  eddy  flux  within  a  depth  of  100  meters  is  usually  less  than  ten  per¬ 
cent  of  the  total  flux. 

Based  upon  thiB  empirical  evidence,  it  has  been  found  useful  to  employ  the  approxi¬ 
mation  that  the  eddy  flux  of  heat,  vapor,  and  momentum  is  constant  throughout  the  layer 
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in  contact  with  the  ground  surface.  The  layer  is  referred  to  as  the  “contact  layer”  in 
this  paper. 


10.  Contact  Laver 

The  momentum,  heat  and  vapor  fluxes  are  assumed  to  be  non-divergent  within  a 
layer  of  50  m  thickness  in  contact  with  the  Earth’s  surface.  The  value  of  the  temperature 
and  specific  humidity  are  computed  over  land  at  the  level  of  the  instrument  shelter  follow¬ 
ing  the  method  outlined  in  Sections  19  and  20.  The  wind  speed  is  assumed  to  vanish  at 
the  level  of  the  surface  roughness.  The  stress  acting  at  the  surface  is  estimated  from 
the  surface  Rossby  number  and  the  static  stability  using  the  empirical  analysis  of 
Lettau  [17]  and  Blackadar  [3].  The  mixing  coefficients  are  computed  using  the  appropriate 
formulas  suggested  by  the  theories  of  Monin  [23]  and  Priestley  [25].  The  eddy  Prandtl 
number  is  taken  to  be  a  constant:  =  K^,  in  forced  convection,  and  =  1.3  K^,  in 

free  convection.  When  free  convection  is  indicated  to  prevail  in  the  contact  layer,  we 
continue  to  utilize  the  forced  convection  wind  profile  between  the  surface  roughness 
height  and  1  m.  A  minus  1/3  power  law  applies  throughout  the  remainder  of  the  contact 
layer. 

The  role  of  the  contact  layer  within  the  model  is  to  provide  lower  boundary  condi¬ 
tions  for  the  equations  governing  the  model  atmosphere  in  the  deeper  transition  layer. 

For  temperature  and  specific  humidity,  we  compute  the  eddy  flux  within  the  contact  layer 
and  use  this  quantity  as  a  boundary  condition.  For  the  wind  equations,  we  use  the  wind 
profile  equations  to  specify  the  wind  speed  and  the  angle  of  geostrophic  deviation  at  the 
base  of  the  transition  layer.  We  assume  that  no  eddy  flux  of  liquid  water  occurs  within 
the  contact  layer. 

11.  Surface  Stress 

The  friction  velocity,  u+,  is  related  to  the  surface  stress  t^,  and  the  air  density, 
p,  by  the  expression 


(11-27) 


Lettau  [17]  has  shown  that  if  one  defines  a  geostrophic  drag  coefficient,  C,  by 
C  =  u*/  G 


(11-28) 
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in  which  G  is  the  surface  geostrophic  wind  speed,  then  an  empirical  relationship  between 
C  and  the  surface  Rossby  number,  R^,  may  be  derived  from  observational  data.  R^  is 
defined  in  terms  of  G,  the  Coriolis  parameter,  f,  and  the  surface  roughness,  Z^,  by 
RQ  =  G  /  (fZQ) .  (11-29) 

In  Lettau’s  analysis,  u*  was  evaluated  for  neutral  stratification.  When  the  temperature 
lapse  rate  departs  from  neutral,  systematic  departures  from  the  previous  estimates  of 
C  occur.  For  lapse  stratification,  the  new  drag  coefficient  is  about  twenty  percent 
greater;  for  moderately  strong  inversion  conditions,  the  new  drag  coefficient  is  some 
twenty  percent  smaller. 

Using  the  data  presented  in  Blackadar’s  paper  [3],  we  obtained,  by  least-squares 
fitting,  the  relation, 

u*  =  G(0. 07625  -  0.00625  log  R  ) .  (11-30) 

n  0 

The  subscript,  n,  indicates  that  this  estimate  is  appropriate  for  neutral  conditions. 
According  to  Blackadar,  some  deviation  of  u+  about  the  value  given  in  the  equation  above 
can  be  attributed  to  the  presence  of  a  geostrophic  wind-shear  within  the  boundary  layer. 
The  fitted  relation  [Eq.  (11-30)]  agrees  very  closely  with  Lettau’s  analysis  [17,  Fig.  3, 

P,  246]. 

12.  Geostrophic  Deviation  Angle 

The  angle  of  deviation  between  the  surface  geostrophic  wind  and  the  wind  in  the 
contact  layer  was  shown  to  be  a  function  of  the  surface  Rossby  number,  R^,  by  Blackadar 
[3].  Using  the  method  of  least  squares  we  fitted  the  following  expression  to  his  data; 

=  a  [log  Rq]2  +  b  log  Rq  +  c  (U-31) 

where  ip  is  the  deviation  angle  in  degrees  and  the  logarithm  is  to  the  base  10.  The 
coefficients  were  computed  to  be, 
a  =  0.625 
b  =  -12.750 
c  =  80.625 

5  10 

ip  varies  between  32.5"  for  R^  =  10'  and  15.6°  for  R^  =  10 
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13.  Turbulence  Regimes 

We  have  used  the  expressions  forced  and  free  convection.  In  forced  convection, 
both  inertial  and  buoyancy  forces  act  to  promote  turbulent  eddying  of  the  air.  In  free 
convection,  buoyancy  forces  are  alone  significant. 

At  all  times  in  the  real  atmosphere  there  is  some  air  motion  and,  as  a  consequence 
of  viscosity,  there  must  always  be  significant  inertial  influence  present  close  to  the 
air— ground  interface.  With  strong  insolation,  it  is  quite  common  under  inactive  synoptic 
conditions  to  find  very  strong  thermal  lapse  conditions  near  the  ground.  In  these  cases, 
the  role  of  the  buoyancy  forces  may  dominate  in  the  excitation  and  maintenance  of 
turbulent  heat  exchange,  but  only  at  some  distance  above  the  interface. 

We  have  found  it  to  be  necessary  to  allow  for  the  occurrence  of  the  free  convective 
regime  within  the  model.  Within  the  contact  layer,  the  transition  to  free  convection  is 
assumed  to  occur  throughout  its  depth  when  the  Richardson’s  number  reaches  a  value 
near  -0.03  [25]. 


14.  Forced  Convection  Relations 

Under  forced  convective  conditions,  Monin’s  similarity  theory  [23]  may  be  shown 
[24]  to  be  consistent  with  the  following  expression  for  the  mixing  coefficient,  K, 

2  8S 

K  =  [kZ  (1  -/JR.]  —  (11-32) 

where 

k  is  von  Karman’s  constant 
Z  is  height 

/3  is  an  empirical  constant 
R.  is  Richardson’s  number 

l 

S  is  the  wind  speed 

In  the  derivation  of  this  expression,  the  mixing  coefficient  was  assumed  to  be  identical 
for  both  heat  and  momentum. 

The  assumption  of  constant  flux  may  be  expressed  in  the  equations, 


-  _H_  „  KS§ 

V  92 


=  U*9* 


(11-33) 
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(11-34) 


-  p  =  KH  =  u*q 
p  =  Kf =  u* 


> 


where 


H 

C 

P 

P 

0 

u* 

9* 

Q 

q 

q* 


is  the  eddy  heat  flux 

is  specific  heat  at  constant  pressure  of  air 
is  the  air  density 
is  the  potential  temperature 
is  the  friction  velocity 

is  a  constant  with  dimensions  of  temperature 

is  the  vapor  flux 

is  the  specific  humidity 

is  a  dimensionless  constant  (gm/gm) 


Using  the  approximation 
90  =  9T  j 

az  az  +  1 


where  y  is  the  dry-adiabatic  lapse  rate,  and  the  boundary  conditions, 


q  =  q. 

at 

z  =  z., 

i 

1 

T  =  T. 

at 

z  =  z. , 

X 

1 

S  =  0 

at 

z  =  z0, 

we  may  integrate  the  constant  flux  equations  to  obtain 


u.  I  z  I  9  A 

S(Z)  =  rto  r  +  <z-zo>’ 


* 

0, 


s,  .  z. 

T(Z)  =  T[  -  r  (Z-Z.)  +  T  to  F U  —  (Z-Z,) 

'  i'  Su* 


q(Z)  =  q.  +  4-  fn  I  + 


0gq*e* 

Su-l 


(Z-Z.) 


(11-35) 


(11-36) 


(11-37) 


(11-38) 

(11-39) 

(11-40) 


in  which  0  is  the  mean  temperature  of  the  layer,  which  is  taken  to  be  a  constant. 


12 


If,  at  the  top  of  the  contact  layer,  the  values  of  T  and  q,  say  T,  and  q,  ,  are 

h  h 

known,  we  may  solve  the  profile  equations  for  0+  and  q+.  Using  the  temperature 
profile,  one  has  a  quadratic  equation  in  0+.  Of  the  two  roots,  the  correct  value  for 
neutral  stratification  is  given  by  only  one  root, 


-u*  S  In  (h/Z.) 
2/3gk(h  -  Z.) 


"k(h  -  Z.f 

2T/21 

£n  (h/Z.) 

)  J 

From  the  profile  equation  for  q,  one  easily  obtains  the  solution  for  q+ 


q*  =  (qh  "  qi) 


In  [h/Z  ]  + 

- ; -  +  [h  -  Z.l 

k  u*  0  1 


(11-41) 


(11-42) 


The  wind  speed  at  the  top  of  the  contact  layer,  S^,  may  be  readily  evaluated  from 


the  wind  speed  profile  equation. 


U*  , 

S.  =  ~fn  — 


e*/3g 


Su„ 


<h  "  Zo> 


(11-43) 


15.  Free  Convection  Formulas 

Priestley’s  similarity  theory  for  free  convection  [25]  yields  the  following  expres¬ 
sion  for  the  mixing  coefficient  for  heat, 


in  which  X  is  an  empirical  constant. 

Utilizing  this  expression  in  the  equation  for  constant  heat  flux,  one  obtains, 


-  X  Z 


£ 

s, 


1/2 


9T 

az 


3/2 


=  u*9* 


(H-45) 


This  expression  may  be  integrated  using  the  boundary  condition,  T  =  T  at  Z  =  Z^, 
to  get  the  profile  formula 


T(Z)  =  T.  -  y  (Z-Z.) 


3u*9* 


x2/3lu,e,|1/3 


(11-46) 
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If  one  assumes  that  the  mixing  coefficient  for  water  vapor  is  identical  to  that  for  heat, 
the  specific  humidity  profile  equation  becomes, 


3u*<l* 


q(Z)  =  qi  '  “W 


Ke*i 


1/3 


(11-47) 


Under  the  assumption  that  the  mixing  coefficient  for  heat  is  thirty  percent  larger 
than  that  for  momentum,  one  may  derive  a  wind  profile  for  free  convection,  vis., 


3.9u, 


S(Z)  =  - 


1 3*1 1/3  L "1/3  „  -1/3 


x2/3  |u*e 


* * 


,1/3  U I  U 


-  z. 


l] 


(11-48) 


We  have  adjusted  this  equation  by  assuming  that  a  near-neutral  forced  convective 
regime  exists  between  and  Z^  +  1  m.  This  yields  the  wind  profile, 


**  |  Z  +  lm 

S(Z)=  Tc^(-V~ 


3.9u„ 


0|1/3 


X2/3|u,0J1/3'g 


(11-49) 


F 


Z"1/3  -  (ZQ  +  1  m)"1/3 


1 


Now  if  the  temperature  and  specific  humidity  are  known  at  Z  =  h,  the  top  of  the 
contact  layer,  we  may  solve  the  profile  equations  for  9+  and  q+, 


The  wind  speed  at  the  top  of  the  contact  layer  is  readily  computed  to  be  given  by 


U* 

fV1”!  3-9u*  | 

lei1/3 

S(h)  : 

=  — -  In 
k 

i  Zo  /  x2/3|Ul(!ej1/3 

li  1 

(11-52) 

CO 

> 

.  V  . 

- 
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1G.  Determination  of  Regime 

The  two  sets  of  results  derived  above  may  be  applied,  provided  that  the  appropriate 
convective  regime  is  determined.  We  discovered  that  the  forced  convection  formula  for 
0+,  which  involves  a  radical,  can  be  used  to  determine  the  character  of  the  regime.  We 
note  that  the  quantity  G+  will  be  complex  if 


h  -  Z. 


r  < 


uif 

4/3g 


fn(  h/Z.) 


k(h 


zi> 


(11-53) 


Upon  utilizing  the  forced  convection  heat  flux,  evaluated  at  a  temperature  difference 
for  which  the  radical  vanishes  in  the  formula  for  9+,  we  can  evaluate  the  Richardson’s 
number  (R.)  as  a  function  of  height,  and  the  parameter,  /3.  Using  h  =  50  m,  Z.  =  1  m, 

/3  =  2.0,  and  Z  =  1.5  m,  we  found  [12]  that  R^  =  -0.032.  This  value  agrees  well  with  the 
critical  Richardson’s  number  quoted  by  Priestley  [25].  We  concluded  that  the  satisfaction 
of  the  inequality  above  was  a  sufficient  condition  for  deciding  that  a  free  convection  regime 
prevailed. 

A  further  test  of  regime  was  subsequently  developed  from  the  formula  for  the  forced 
convection  mixing  coefficient  evaluated  at  the  top  of  the  contact  layer, 

0g1 


■S. 


=  hku. 


1.0  + 


;khu*e*“l 

J 


(11-54) 


This  is  equivalent  to 

=  ku*h  (1  -  f3  R.)  (Z  =  h) 


(11-55) 


and  it  is  clear  that  the  factor  (1  -  ^R^)  must  be  positive.  This  leads  to  a  condition 


P gkhu^e* 

1.0  +  - - —  >  0 

Su" 


(11-56) 


The  existence  of  this  added  condition  reflects  the  possible  existence  of  a  mixture 
of  free  and  forced  convective  regimes  within  the  contact  layer  of  a  more  complex  nature 
than  that  indicated  in  our  treatment  of  the  wind  profile  in  free  convection.  The  description 
of  this  situation  is  beyond  the  scope  of  our  model.  We  have,  of  course,  been  forced  to 
recognize  the  constraint  imposed  by  the  last  inequality,  which  is  more  stringent  than  the 
inequality  given  earlier  In  practice,  we  have  gone  a  step  further  and  employed  the  free 
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convection  formulas  in  lapse  conditions  whenever  the  free  convection  heat  flux  exceeds 
that  given  by  the  forced  convection  formula. 

17.  Continuity  of  Heat  Flux 

It  seemed  reasonable  to  require  the  heat  flux  to  be  continuous  at  the  critical 
Richardson’s  number.  We  found  [12]  that  this  condition  was  adequate  to  relate  the  value 
of  the  parameter,  X,  to  the  value  of  p.  The  resulting  expression  is  of  the  form 

X  =  k2  c  (11-57) 

The  symbol,  c,  stands  for  a  constant  dependent  solely  upon  the  levels  Z.  and  h. 

Using  k  =  0.38,  we  found 

X  =  0.85  v?  (11-58) 

For  p  =  2.0,  X  =  1.2,  a  value  between  those  given  for  this  parameter  by  Priestley  [25] 
and  Dyer  [6]. 


18.  Extreme  Stability 

In  addition  to  the  limit  on  the  applicability  of  the  Monin  forced  convection  formulas 

in  lapse  conditions  indicated  above,  we  note  that  the  formulas  become  invalid  when  the 

Richardson  number  becomes  larger  than  p  1.  In  this  case,  we  assume  that  the  quantities 

0*  and  qt  may  be  computed  using  the  constant  minimal  value  for  the  mixing  coefficient, 

2  2-1 

K  .  (~10  cm  sec  ). 

mm 


In  this  way,  one  obtains, 


K 


= 


min 

u* 


and 


Tl  -  T. 
h  i 

h  -  Z. 


+  y 


(11-59) 


K  .  fq.  -  q.“| 

_  mm  I  h  l  I 

u.  Lh  -  Z.J 


(11-60) 


To  compute  the  wind  speed  at  Z  =  h,  we  assume  that  the  wind  is  a  fixed  fraction  of  the 
surface  geostrophic  wind  speed,  G, 


S(h)  =  0.176  G  .  (11-61) 

The  constant  (0.176)  was  estimated  from  the  “Ekman  Spiral”  using  K  =  104  cm^ 
sec  1,  f  =  10"4  sec-1,  and  the  assumption  that  the  wind  vanishes  at  Z  =  0  and  becomes 
geostrophic  as  Z  —  00 . 
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19.  Surface  Temperature  Specification 

By  surface  temperature,  we  mean  air  temperature  measured  within  an  instrument 
shelter  located  1.22  m  above  the  ground.  In  the  case  of  points  located  over  the  ocean, 
we  consider  the  sea  surface  temperature  to  be  the  appropriate  surface  temperature. 

The  temporal  variation  of  the  surface  temperature  over  the  land  must  be  computed 
within  the  model.  This  variation  can  be  attributed  to  four  physical  processes:  divergence 
of  radiative  heat  flux,  divergence  of  eddy  heat  flux,  thermal  advection,  and  latent  heat 
exchange.  In  an  effort  to  evaluate  the  influence  of  the  first  two  processes,  we  made  use 
of  an  empirical  technique  for  specifying  the  temporal  variation  of  surface  temperature. 
This  technique,  developed  by  Bryan  [5],  involves  the  use  of  the  equation 


in  which 


^-T(t)  +  bxT(t)  =  bQ  +  b2  S(t)  +  b3r(t) 


S(t) 

S(t) 

r(t) 

r(t) 


7Tt 

sin  6  sin  <p  -  cos  6  cos  (p  cos  12  (R  s  t  SS) 

0  (otherwise) 

.  .  7Tt 

(9/9t)  S(t)  =  X/12  cos  <5  cos  <p  sin  — 

JL  ct 

(R  2£  t  £  12) 


0  (otherwise) 


(11-62) 


(11-63) 


(11-64) 


and  <5  is  the  solar  declination,  <p  is  the  latitude,  R  is  the  local  time  of  sunrise  in  hours, 

S  is  the  local  time  of  sunset  in  hours,  T  is  the  temperature,  and  t  is  the  time  in  hours 
after  local  midnight.  The  coefficients,  b  ,  b  ,  b  ,  and  b  were  estimated  using  station 

U  X  Z  o 

records  of  hourly  temperature  change.  The  coefficients  were  categorized  by  month  and 
by  cloud  cover.  The  basic  temperature-change  data  were  derived  using  ten -year  average 
diurnal  temperatures  from  which  any  net  diurnal  tendency  (attributable  to  advection)  was 
first  subtracted.  The  data  used  in  our  computations  of  the  coefficients  were  obtained 
from  the  U.S.  Air  Force  and  are  described  in  a  technical  note  by  Kimball,  Richardson, 
and  Frey  [15]. 

Once  the  empirical  coefficients  have  been  derived,  the  tendency  implied  by  Bryan’s 
equation  may  be  computed  given  the  type  of  cloudiness,  the  present  temperature,  the 
local  time,  geographical  position  and  the  time  of  year.  The  cloudiness  specification  was 
chosen  to  take  advantage  of  the  fact  that  low  cloudiness  would  be  predicted  within  the 
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model.  Middle  or  high  cloudiness  might  reasonably  be  predicted  by  a  free  air  model 
and  made  available  to  this  model.  The  three  categories  of  cloudiness  were  therefore 
taken  to  be 

(a)  Category  1— clear  or  scattered  clouds 

(b)  Category  2— broken  or  overcast  clouds  above  5000  ft 

(c)  Category  3— broken  or  overcast  clouds  below  5000  ft. 

Thermal  advection  at  the  level  of  the  instrument  shelter  cannot  be  neglected.  We 
compute  it  by  assuming  a  logarithmic  wind  profile  to  exist  through  the  1.22-m  air  layer 
above  the  surface  roughness  height. 

The  influence  of  latent  heat  exchange  is  neglected  as  a  smaller  order  effect  in  this 
model. 

20.  Surface  Humidity  Specification 

The  model  requires  the  specification  of  the  temporal  variation  of  specific  humidity 
at  the  lower  boundary  of  the  contact  layer.  For  points  which  lie  over  water,  we  may  safely 
assume  that  the  relative  humidity  is  always  close  to  one  hundred  percent  near  the  air- 
water  interface.  Over  land,  on  the  other  hand,  the  low-level  relative  humidity  depends 
upon  the  amount  of  available  soil  moisture  and  the  complex  process  by  which  this  moisture 
is  liberated  from  the  soil. 

We  originally  considered  approaching  the  prediction  of  surface  relative  humidity 
in  a  manner  analogous  to  that  used  for  the  surface  temperature.  The  Air  Force  data 
collection  provides  the  requisite  dew-point  temperatures  for  such  an  approach.  It  was 
pointed  out  to  us  that  another  technique  was  available  which  would  not  require  any 
significant  development  effort  and  we  therefore  decided  to  adopt  it  for  the  present. 

In  their  effort  to  develop  an  analog  computer  for  micrometeorological  use, 

Halstead,  et  al.  [13]  were  faced  with  the  need  to  compute  the  percentage  of  available 
energy  used  to  evaporate  water.  They  employed  a  parameter,  M,  which  is  denoted  as 
the  “percent  wetted  area”  in  physiological  climatology.  For  their  problem,  Halstead, 
et  al.,  related  M  to  the  water  vapor  densities  as  follows, 

p'o  -  p'h  ■  M  <p'oSAT  -  V  <u-65> 

where  p'  denotes  water  vapor  density;  the  subscript,  0,  denotes  a  measurement  near 
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the  surface,  the  subscript,  H,  denotes  a  measurement  some  distance  removed  from  the 

surface;  and  p'q  ,  is  the  saturation  vapor  density  near  the  surface.  Using  mete- 
SAT 

orological  data  gathered  during  the  Great  Plains  field  program  (this  program  is 
described  by  Lettau  and  Davidson  f  18] X  it  was  found  [13]  that  M  was  well  correlated  to 
measured  soil  moisture  and  that  its  value  tended  to  remain  constant  during  periods  of 
twenty-four  hours  or  more. 

We  have  adapted  this  result  to  our  problem  as  follows:  let  q  be  the  specific 

humidity  at  the  top  of  the  contact  layer,  and  be  the  surface  specific  humidity  measured 

at  the  level  of  the  instrument  shelter;  finally,  let  q.  be  the  saturation  specific  humidity 

1  s 

measured  at  the  same  level  as  q..  Using  measurements  of  these  quantities  made  prior 
to  the  initial  time  of  the  forecast  we  may  compute  a  value  of  M  from 


M  =  (q  -  q  )  /  (q  -  q  ).  (11-66) 

1  IB 

Now,  holding  M  constant  through  the  forecast  interval  permits  one  to  solve  for  the 
surface  relative  humidity  from  the  equation, 


RHS(t)  =  M  +  (1  -  M)  q(t)  /  q.  (t) .  (11-67) 

IS 

Since  q.  is  principally  a  function  of  temperature,  its  value  may  be  computed  from  the 

IS 

predicted  value  of  the  surface  temperature.  The  value  of  q  is  predicted  within  the 
model  and  is  available  at  every  time  step  for  use  in  evaluating  Eq.  (11-67)  for  RHS(t). 


21.  Coordinate  System  and  Terrain  Variation 

The  equations  used  in  this  prediction  model  are  written  in  Cartesian  coordinates. 
The  horizontal  coordinates  (x,  y)  are  established  on  a  polar  stereographic  mapping  of 
the  Northern  Hemisphere  with  their  origin  taken  at  the  North  Pole.  The  particular  polar 
stcreographic  map  used  is  “true”  at  60°  North  Latitude  and  has  a  scale  of  1:5,000,000. 
The  map  scale  factor,  a,  is  defined  as  follows, 

u  =  1 1  +  sin  ~  ]  /  [  1  +  sin  <p]  (11-68) 

O 

and  its  numerical  value,  over  the  region  shown  in  Fig.  1,  is  given  in  Table  I. 

The  vertical  coordinate  (Z)  has  its  origin  at  the  height  of  the  surface  terrain.  The 
elevation  of  the  surface  terrain  was  derived  from  the  values  given  by  Berkofsky  and 
Bertoni  [4]  and  from  other  sources  of  topographic  data.  The  elevations  assigned  to  the 
the  grid  points  are  given  in  Table  II. 
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Fig.  1.  Grid  network  used  in  this  study. 
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TABLE  I 

MAP  SCALE  FACTOR  AT  GRID  POINTS  IN  FIG.  1 


\l 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

10 

1.104 

1.100 

1.099 

1.097 

1.095 

1.095 

1.094 

1.095 

1.095 

1.096 

9 

1.115 

1.114 

1.111 

1.109 

1.108 

1.108 

1.108 

1.108 

1.108 

1.109 

8 

1.129 

1.126 

1.124 

1.122 

1.121 

1.120 

1.120 

1.120 

1.121 

1.123 

7 

1.142 

1.140 

1.138 

1.137 

1.135 

1.134 

1.134 

1.135 

1.135 

1.137 

6 

1.156 

1.154 

1.152 

1.151 

1.149 

1.148 

1.147 

1.148 

1.149 

1.151 

5 

1.172 

1.168 

1.167 

1.165 

1.164 

1.163 

1.163 

1.163 

1.164 

1.165 

4 

1.187 

1.184 

1.183 

1.180 

1.179 

1.178 

1.178 

1.178 

1.179 

1.180 

3 

1.201 

1.199 

1.197 

1.196 

1.194 

1.193 

1.194 

1.194 

1.194 

1.196 

2 

1.217 

1.215 

1.213 

1.212 

1.210 

1.209 

1.209 

1.209 

1.209 

1.210 

1 

1.234 

1.232 

1.230 

1.229 

1.228 

1.227 

1.226 

1.227 

1.228 

1.228 

TABLE  B 

ELEVATION  OF  TERRAIN  (IN  METERS)  ABOVE  MEAN  SEA  LEVEL  AT 
VARIOUS  GRID  POINTS  GIVEN  IN  FIG.  1 


L 

M  N. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

10 

367 

206 

227 

191 

227 

191 

191 

185 

70 

82 

9 

267 

327 

212 

191 

245 

191 

182 

167 

121 

182 

8 

261 

203 

197 

233 

258 

191 

348 

667 

369 

124 

7 

221 

179 

185 

245 

306 

282 

370 

385 

82 

15 

6 

252 

109 

118 

158 

273 

173 

606 

148 

9 

0 

5 

124 

91 

148 

239 

336 

439 

364 

124 

0 

0 

4 

70 

115 

183 

211 

330 

530 

221 

70 

6 

0 

3 

48 

70 

185 

203 

324 

136 

58 

15 

0 

0 

2 

115 

92 

52 

140 

93 

39 

18 

0 

0 

0 

1 

61 

52 

39 

61 

70 

21 

0 

0 

0 

0 

The  surface  roughness  parameter  has  been  derived  from  data  given  by  Kung  [16], 
He  made  use  of  information  on  land  usage  available  from  the  U.S.  Department  of 
Agriculture  and  related  it  to  estimates  of  the  surface  roughness  parameter,  available 
for  various  crops  and  trees,  from  detailed  local  studies.  The  grid-point  values  of  sur¬ 
face  roughness  used  in  the  model  integrations  of  winter  cases  are  given  in  Table  HI. 
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TABLE  HI 

SURFACE  ROUGHNESS  PARAMETER  (IN  CENTIMETERS)  AT  THE 
GRID  POINTS  SHOWN  IN  FIG.  1 


\ 

M  \ 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

10 

1 

1 

1 

1 

2 

1 

2 

10 

12 

32 

9 

1 

1 

1 

1 

2 

1 

1 

2 

22  * 

32 

8 

5 

2 

2 

2 

2 

1 

2 

19 

30 

20 

7 

10 

5 

4 

3 

2 

5 

15 

32 

21 

10 

6 

15 

15 

15 

15 

16 

20 

30 

32 

10 

1 

5 

25 

28 

30 

30 

30 

35 

45 

30 

10 

1 

4 

35 

40 

42 

50 

52 

61 

47 

20 

10 

1 

3 

45 

50 

56 

60 

70 

65 

25 

10 

1 

1 

2 

55 

62 

70 

80 

93 

50 

8 

1 

1 

1 

1 

70 

90 

100 

101 

100 

20 

1 

1 

1 

1 
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SECTION  III 


INPUT  DATA  ANALYSIS 

The  surface  and  upper-air  observations  made  during  two  time  periods  were  col¬ 
lected  through  local  weather  teletype  facilities  with  the  cooperation  of  the  personnel  of 
the  Travelers  Weather  Service.  From  these  routine  data,  we  constructed  the  grid-point 
data  required  by  the  prediction  model  We  employed  some  machine  computation  methods 
in  the  data  preparation  but  the  analysis  procedure  was  not  fully  automated.  The  surface 
charts  were  analyzed  for  the  temperature  field.  Grid  point  values  of  surface  tempera¬ 
ture  were  interpolated  from  the  surface  isopleths.  Over  the  ocean,  we  used  monthly 
mean  values  of  the  sea  surface  temperature  made  available  to  us  by  Mr.  A.  Thomasell 
of  the  Research  Center  staff. 

Radiosonde  reports  for  each  observing  station  were  processed  to  obtain  tempera¬ 
ture  and  dew-point  temperature  values  at  each  level  in  the  vertical  required  in  the  model. 
These  values  were  linearly  interpolated  from  mandatory-  and  significant-level  observations. 
We  then  plotted,  at  the  various  reporting  stations,  the  difference  in  temperature  between 
particular  levels.  These  data  were  then  analyzed  and,  by  interpolation  between  isopleths, 
grid-point  values  were  obtained  When  these  vertical  temperature  differences  were  used 
in  conjunction  with  the  grid-point  values  of  surface  temperature,  we  reconstructed  vertical 
temperature  distributions  at  each  grid  point. 

From  the  interpolated  radiosonde  data,  we  plotted  the  dew-point  depression  at  each 
reporting  station.  Analyses  of  these  fields  were  then  constructed  and  interpolation  again 
yielded  grid-point  values.  When  these  data  were  combined  with  the  grid-point  tempera¬ 
tures,  the  grid-point  values  of  dew-point  temperature  resulted.  Conversion  of  these  data 
to  specific  humidity  values  was  carried  out  using  a  standard  atmosphere  pressure- 
height  relation. 

This  analysis  technique  seems  suitably  controlled,  but  it  is  necessary  to  use  good 
judgement  when  performing  the  hand  analysis.  The  development  of  a  good  automated 
procedure  is  highly  desirable. 

The  geostrophic  wind  components  at  Z  =  II  and  the  upper  level  cloudiness  can,  in 
principle,  be  obtained  from  free-air  prediction  models.  However,  we  were  forced  to  use 
the  observed  synoptic  reports  to  construct  the  needed  input. 
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Geostrophic  winds  were  computed  from  the  observed  field  of  geopotential  at  850 
and  700  mb  at  twelve  hour  intervals.  These  were  interpolated  vertically  to  obtain  geo¬ 
strophic  winds  at  Z  =  H.  Six-hour  pibal  wind  observations  were  used  to  determine 
if  the  temporal  variation  of  the  geostrophic  winds  between  the  twelve-hour  observations 
was  reasonably  accurate.  Changes  were  introduced  at  the  sixth  hour  if  it  appeared 
necessary. 

Nephanalysis  of  the  middle-  and  high-cloud  fields  was  attempted  using  the  six- 
hourly  surface  synoptic  data.  The  mean  condition  of  cloudiness  was  fixed  over  a  six- 
hour  interval  at  each  grid  point  by  assigning  to  the  parameters  ICLU1  and  ICLU2  the 
value  1  for  clear  or  scattered  high  and  middle  cloudiness  or  2  for  broken  to  overcast 
high  or  middle  cloudiness. 

It  should  be  noted  that  we  did  not  attempt  to  modify  the  radiosonde  temperature 
or  humidity  observations  to  reflect  the  observed  cloudiness.  Such  a  procedure  might 
prove  desirable  in  operational  use  of  the  model.  In  this  connection  the  initial  specific 
moisture  (r)  field  was  always  equal  to  the  initial  specific  humidity  (q)  field. 
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SECTION  IV 


ANALYSIS  OF  TESTS 

The  results  of  three  forecasts  (identified  as  Cases  A,  B,  and  C)  made  with  the 
prediction  model  are  given  below.  Each  test  case  is  a  12-hr  forecast. 

Case  A:  12Z  February  6,  1964  to  OOZ  February  7,  1964 

This  case  was  run  previously  with  an  earlier  model  [12]  and  served  as  a  control  in 
the  development  of  the  revised  model.  The  synoptic  charts  at  the  initial  time  [see  Fig.  2(a) 
to  2(n)]  depicted  a  diffuse  occluded  low  at  the  surface  with  its  center  in  southern  Ohio.  At 
upper  levels,  this  low  was  more  intense  but  the  cold  air  was  located  to  its  south.  On  the 
surface  chart  at  the  initial  time  “four-dot  rain”  was  reported  at  Cape  Hatteras,  N.  C. 
and  Norfolk,  Va.,  indicative  of  the  development  of  a  secondary  cyclone.  This  secondary 
circulation  was  very  evident  in  the  50-m  wind  field  diagnosed  using  the  model  (see  Fig.  3). 
During  the  forecast  period,  the  secondary  circulation  intensified  and  moved  rapidly  north¬ 
ward.  At  upper  levels  this  development  was  associated  (during  the  first  six  hours)  with  a 
gradual  turn  of  the  upper  winds  from  southerly  to  westerly  over  the  south-central  portion 
of  the  region.  During  the  second  half  of  the  forecast  period  the  upper  circulation  rapidly 
adjusted  to  the  developing  secondary  cyclone.  The  upper  low  accelerated,  and  at  the  end 
of  the  forecast  period  was  located  over  eastern  Pennsylvania. 

The  results  obtained  in  this  case  [see  Fig.  2(a)  through  2(n)]  can  be  seen  to  be  of 
highest  quality  in  regions  removed  from  the  lateral  boundaries.  This  is  to  be  expected  in 
view  of  our  neglect  of  horizontal  advection  on  inflow  boundaries.  The  most  serious  errors 
introduced  in  this  case  occur  in  the  northeastern  region.  At  the  initial  time,  the  air  in 
this  region  was  very  dry.  As  the  situation  evolved,  the  humidity  increased  due  to  the  over¬ 
water  trajectory  of  the  air  in  advance  of  the  secondary  cyclone.  The  model  does  not  detect 
this  transport  because  of  the  boundary  condition.  The  movement  of  the  cold  air  trough 
through  the  center  of  the  region  is  well  predicted,  as  is  the  cloudiness  and  precipitation 
observed  through  the  Appalachians.  The  close  correspondence  of  the  observed  and  pre¬ 
dicted  50-m  wind  fields  bears  witness  to  the  adequacy  of  the  predicted  thermal  field  con¬ 
figuration.  Since  the  thermal  wind  within  the  boundary  layer  is  predicted,  such  close 
agreement  of  the  50-m  wind  fields  could  not  occur  if  the  predicted  thermal  field  was  poorly 
correlated  with  the  observed  field. 
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The  pattern  of  humidity  (cloudiness)  depicted  in  the  vertical  cross-sections  is 
traceable  to  the  ‘‘terrain  induced”  component  of  the  vertical  wind.  In  Fig.  4,  we  show 
this  field  as  predicted  at  the  verification  time.  It  is  significant  that  the  scale  of  the 
major  feature  of  this  field  is  smaller  than  the  NMC  mesh  spacing.  The  computation  of 
the  terrain  influence  on  the  weather  associated  with  the  vertical  velocity  field  requires 
the  consideration  of  such  smaller-scale  phenomena.  Indeed,  the  secondary  cyclone  at  the 
initial  time  (see  Fig.  3)  also  would  have  been  lost  in  the  NMC  grid. 

Figures  (5)  and  (6)  show,  respectively,  the  net  sensible  heat  and  latent  heat  added 
to  the  transition  layer  during  the  forecast  periods.  The  values  plotted  were  computed  by 
integrating  the  “instantaneous  values”  of  the  fluxes  at  hourly  intervals.  It  is  interesting 
to  note  that  the  maximum  value  of  sensible  heat  transfer  occurs  over  the  land.  Clear  skies 
prevailed  over  this  region  and  the  pattern  of  heating  reflects  the  cloudiness  distribution 
to  a  significant  extent.  Over  the  ocean,  the  latent  heat  transfer  is  much  greater  than  the 
sensible  heat  transfer.  The  accuracy  of  this  result  is  in  doubt  due  to  the  inadequacy  of 
the  initial  data  in  that  region. 

Case  B:  00Z  January  23— 12Z  January  23,  1965 

Figure  7(a)  and  (b)  diqjhysthe  synoptic  patterns  at  the  initial  and  final  times  of  the 
forecast  period.  Low  cloudiness  and  precipitation  are  noted  across  the  northern  and 
western  portions  of  the  region  at  the  initial  time.  The  activity  in  the  north  is  associated 
with  a  stationary  front.  In  the  southwestern  region,  showery  conditions  are  reported  in 
the  warm  air— the  NMC  analysis  indicated  the  presence  of  a  squall  line  in  this  area. 

During  the  forecast  interval,  the  cold  air  pushed  southward  along  the  eastern  boundary, 
with  a  small  depression  indicated  on  the  front  near  Washington,  D.C.  Precipitation  continued 
to  occur  over  the  northern  portion.  Low  cloudiness  developed  over  the  central  portion  of 
the  map  and  an  elongated  area  of  precipitation  was  reported  along  the  western  slope  of  the 
Appalachian  Mountains. 

Figure  7(c),  (d),  and  (e)  shows  the  analysis  of  the  thermal  and  humidity  fields  observed 
at  12Z  January  23,  1965  The  corresponding  forecast  charts  are  given  in  Fig.  7(f),  (g),  and 
(h).  The  isotherm  ribbon  associated  with  the  stationary  front  was  observed  to  move  south¬ 
ward  along  the  east  coast  during  the  forecast  interval.  This  movement  was  not  predicted 
by  the  model  and  the  intensity  of  the  thermal  gradient  in  that  region  is  under-predicted. 
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The  cold  trough  through  the  central  region  is  a  reflection  of  the  terrain  elevation.  Rather 
little  temperature  change  occurred  therein.  The  incursion  of  higher  humidity  into  the 
southwestern  portion  is  predicted  by  the  model,  principally  at  the  upper-most  level  and 
at  the  level  of  the  instrument  shelter.  The  observed  and  predicted  surface  charts 
(instrument  shelter  temperature  and  50-m  wind)  shown  in  Fig.  7(i)  and  (j)  are  in  good 
agreement  except  in  the  northeastern  portion.  The  band  of  precipitation  on  the  western 
slope  of  the  Appalachian  Mountains  seems  to  be  due  to  a  high-level  condensation  process. 
The  low  cloud  reported  is  apparently  due  to  the  advective  formation  of  stratus  and, 
possibly,  in  part  due  to  the  evaporation  of  precipitation  falling  into  the  boundary  layer 
from  above. 

West  to  east,  vertical  cross  sections  through  the  predicted  atmosphere  are  given 
in  Fig.  7(k)  through  (n).  In  Fig.  8  we  have  constructed  a  south  to  north  vertical  cross 
section  along  the  coordinate  line,  L  =  5. 

Figures  9  through  13  are  soundings  showing  the  initial  and  predicted  grid-point 
temperature  structure  and  the  thermal  structure  at  the  verification  time  observed  at  a 
nearby  radiosonde  station.  In  Fig.  9  we  show  a  point  near  the  eastern  boundary  of  the 
grid.  The  front  was  observed  to  pass  south  of  this  point  during  the  forecast  period.  The 
predicted  sounding  is  too  warm  at  low  levels,  reflecting  the  failure  of  the  model  to 
properly  displace  the  frontal  position.  Above  300  m  the  predicted  sounding  agrees  closely 
with  the  observed  data.  We  may  note  also  the  very  strong  inversion  in  the  lowest  50  m  of 
the  predicted  sounding.  This  low  level  structure  may  be  unrealistic  due  to  the  absence  of 
sufficient  refinement  in  the  contact  layer  theory  for  strong  inversion  conditions. 

Figure  10  shows  a  grid  point  near  the  elongated  axis  of  precipitation  on  the  western 
slope  of  the  Appalachians.  The  negligible  change  observed  in  the  vertical  temperature 
structure  is  properly  treated  by  the  model.  Figure  11  shows  a  point  located  just  to  the 
north  of  the  stationary  front  throughout  the  forecast  period.  At  low  levels  (below  600  m) 
some  cooling  was  observed  to  occur,  whereas  the  model  predicted  some  warming.  Note 
again  the  strong  inversion  in  the  predicted  sounding  between  the  surface  and  50  m. 

Figure  12  is  for  a  point  located  similarly  to  that  shown  in  Fig.  11,  with  rather  similar 
results.  In  Fig  13  we  indicate  the  results  at  a  point  located  to  the  north  of  the  front  but 
nearer  to  the  developing  storm  (just  off  the  western  edge  of  the  map).  The  warming  above 
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the  frontal  inversion  is  predicted  by  the  model,  but  the  observed  low-level  cooling  is  not 
indicated  in  the  forecast.  We  may  note  the  northward  displacement  of  the  thermal  ribbon 
in  the  horizontal  depiction  charts  as  the  explanation  for  the  low-level  error  at  this  grid- 
point. 

Figure  14  shows  the  temperature  and  dew-point  soundings  observed  at  radiosonde 
station  226  at  the  initial  and  final  times  of  the  forecast,  as  well  as  the  predicted  values  at 
the  nearby  grid  point.  This  station  was  reporting  moderate  rain  at  the  verification  time. 
The  analysis  is  complicated  by  the  nearness  of  the  grid  point  to  an  inflow  boundary,  but 
with  the  exception  of  the  observed  inversion  between  850  m  and  1150  m  and  the  departure 
below  300  m,  the  predicted  temperature  lapse  rate  is  very  close  to  that  observed.  The 
predicted  dew-point  spread  is  greater  than  the  initial  and  observed  values  above  400  m. 

The  vertical  velocity  predicted  at  this  grid  point  was  small  until  the  tenth  hour.  At  that 
time  an  upward  flow  of  2  cm  sec  *  was  predicted  in  the  upper  portion  of  the  boundary 
layer.  Consequently,  we  must  ascribe  the  predicted  drying  to  advection.  Below  400  m, 
the  observed  increase  of  humidity  is  well  predicted.  In  the  forecast  super-saturation 
was  predicted  between  50  m  and  150  m.  The  liquid  water  concentration  predicted  was 
between  0.2  and  0.4  gms  kg  \  There  is  a  good  possibility  that  the  observed  increase  in 
humidity  at  upper  levels  was  in  response  to  the  evaporation  of  the  precipitation  which 
fell  through  the  boundary  layer  from  above.  One  must  also  note  that  the  convective 
activity  which  occurred  in  this  general  region  is  not  predicted  by  the  physical  model. 

Its  potential  influence  on  the  observed  soundings  cannot  be  ignored  as  a  source  of  spatial 
scale  confusion. 

Case  C:  12Z  January  24,  1965— 00Z  January  25,  1965 

In  Fig.  15(a)  and  (b),  the  synoptic  patterns  at  the  beginning  and  end  of  the  Case  C 
forecast  are  indicated.  A  low-pressure  center  was  located  over  Michigan  at  the  initial 
time.  Along  the  east  coast  a  ridge  of  cold  air  was  entrenched  to  the  east  of  the  Appalachian 
Mountains  and  north  of  the  stationary  front  located  along  the  border  between  Virginia  and 
North  Carolina.  A  cold  front  extended  from  the  low  southward  to  the  southwestern  edge 
of  the  grid  region. 

During  the  twelve  hours  of  the  forecast,  the  low  moved  east— northeastward  to  the 
north  of  Lake  Erie.  The  front  moved  eastward  with  the  low.  It  extended  southward  through 
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western  Pennsylvania  and  Virginia  as  a  warm  type  occlusion.  From  Virginia  southward, 
it  was  a  cold  front  along  which  squalls  were  reported.  A  secondary  low  had  developed  on 
the  stationary  front  over  the  coastal  waters  and  at  00Z  was  located  at  40°N  and  71°W. 

The  entire  grid  region  was  covered  by  broken  or  overcast  low  cloudiness  at  the 
initial  time .  Precipitation  was  widespread  An  extensive  area  of  snow  was  reported  to 
the  west  of  the  low.  An  area  of  rain  showers  was  reported  along  the  cold  front  as  far 
south  as  Nashville,  Term.  Under  the  cold-air  ridge  along  the  east  coast,  freezing  pre¬ 
cipitation  was  reported  over  a  large  area.  Showery  conditions  also  existed  in  the  warm  air 
along  the  coast  and  over  the  coastal  waters. 

At  the  end  of  the  forecast  period,  cloudiness  diminished  over  the  southwestern 
portion  of  the  region  but  continued  in  evidence  along  the  northern  half  and  along  the  east 
coastal  region.  The  precipitation  pattern  changed  in  response  to  the  displacement  of  the 
low  and  the  cold  front.  A  small  region  through  central  New  York  and  Pennsylvania 
separated  the  frontal  precipitation  from  that  in  the  cold  ridge  associated  with  the  off¬ 
shore  low. 

Figure  15(c),  (d),  and  (e)  depicts  the  analysis  of  temperature  and  humidity  observed 
at  00Z,  Jan.  25,  1965  The  analyses  of  the  predicted  fields  are  given  in  Fig.  15(f),  (g),  and 
(h).  At  the  three  levels,  the  forecast  of  temperature  in  the  northeast  is  much  too  high. 

The  explanation  for  this  error  lies  in  the  development  of  the  secondary  coastal  low.  In 
Fig.  16,  the  sounding  at  radiosonde  station  518,  near  grid  point  L  =  10,  M  =  7,  is  shown. 

It  indicates  that  below  1500  m  the  forecast  sounding  is  displaced  toward  higher  tempera¬ 
tures.  The  observed  trend  is  the  reverse.  If  one  compares  the  50-m  wind  fields  shown 
in  Fig.  15(i)  and  (j),  the  influence  of  the  off-shore  development  is  evident.  The  low-level 
wind  field  is  distorted  away  from  verification  because  the  thermal  wind  field  is  inaccu¬ 
rately  predicted. 

The  principal  line  of  cold  air  is  well  positioned  in  the  forecast  charts.  We  must 
note  that  the  verification  analysis  over  the  ocean  is  highly  speculative.  As  a  consequence, 
the  thermal  gradients  in  the  southeast  portion  are  not  reliable.  The  errors  in  the  forecast 
temperature  fields  are  given  in  Table  IV  for  each  grid  point,  excluding  those  over  the 
ocean.  Also  given  are  the  errors  resulting  from  a  persistence  forecast.  The  later 
values  are  the  negative  of  the  observed  temperature  changes. 
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TABLE  IV 

COMPARISON  OF  TEMPERATURE  PREDICTION  ERRORS  (°K)* 
THE  MODEL  VERSUS  PERSISTENCE! 


(a)  At  500  m  level 


(b)  At  1150  m  level 


♦Forecast  temperature  minu6  observed  temperature. 

fModel-forecast  errors  are  entered  in  the  left  of  each  box,  persistence  errors 
are  entered  in  the  right. 
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(c)  At  2000  m  level 
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The  average  absolute  error  in  temperature  at  each  level  is  given  in  Table  V  below. 
We  have  computed  a  value  F  over  the  entire  grid  excluding  only  points  over  the  ocean, 

w 

and  a  value  over  the  interior  points  (L  =  3,  8;  M  =  3,  ....  8)  excluding  points  with¬ 

in  two  intervals  of  a  boundary.  Similar  values  (P^,  P^)  were  computed  for  the  persistence 
forecast. 


TABLE  V 

AVERAGE  ABSOLUTE  ERROR  IN  TEMPERATURE 


Z  (m) 

FW  (°K> 

Fj  <°K) 

pw  <°K> 

^  <°K> 

500 

3.56 

3.22 

4.98 

6.72 

1500 

2.81 

2.17 

5.90 

7.86 

2000 

2.64 

2.19 

4.97 

6.64 

The  verification  of  the  humidity  prediction  is  best  done  by  using  the  observed  low 
cloud  and  precipitation.  The  analyzed  radiosonde  humidity  data  rarely  indicates  in  a 
direct  fashion  the  observed  cloudiness  pattern.  It  should  again  be  noted  that  the  initial 
humidity  data  employed  in  the  forecast  were  not  systematically  modified  to  reflect  the 
initial  observed  cloudiness. 
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In  Fig.  17,  the  regions  reporting  broken  or  overcast  low  cloudiness  are  indicated. 
We  also  show  the  regions  reporting  precipitation  together  with  an  indication  of  the  type 
of  precipitation.  Figure  18  was  constructed  using  the  humidity  predictions  at  the  three 
levels,  500  m,  1150  m,  and  2000  m.  If  we  assume  that  regions  with  supersaturation  at 
one  of  the  three  levels  are  experiencing  precipitation,  then  we  may  note  the  high  correla¬ 
tion  in  the  region  of  frontal  precipitation.  Reference  to  Fig.  15(k)  through  (n),  vertical 
cross  sections  through  the  predicted  atmosphere,  will  indicate  that  Fig.  18  does  not 
depict  the  full  detail  available  in  the  forecast.  By  comparison  with  Fig.  19,  which  is  a 
presentation  of  the  analysis  of  the  observed  radiosonde  humidity  field,  the  composite 
prediction  chart  seems  to  be  superior  as  a  depiction  of  observed  weather  conditions. 

To  provide  further  detail  on  the  skill  of  the  prediction  model  we  have  reproduced 
the  predicted  temperature  soundings  at  five  interior  grid  points  located  near  radiosonde 
stations.  The  charts  (see  Figs.  20  through  24)  indicate  the  initial  and  predicted  tempera¬ 
ture  at  the  grid  point  plus  the  observed  temperature  at  the  nearby  radiosonde  station. 
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Fig.  2(a).  Case  A, observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  initial  time.  Shading  indicates  regions  covered  by  broken  or 
overcast  low  cloudiness  ,  or  experiencing  precipitation  £IM  . 
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Fig.  2(b).  Case  A,  observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  verification  time.  Shading  indicates  regions  covered  by  broken 
or  overcast  low  cloudiness  ,  or  experiencing  precipitation  . 
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Fig.  2(c).  Case  A:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80^  s  rh  <  100% 
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Fig.  2(d).  Case  A:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  so%  s  rh  <  100% 
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Fig.  2(e).  Case  A:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  80%  ^  rh  <  ioo% 
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Fig.  2(f).  Case  A:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80%  s  rh  <  100%  l&gga 

Rii  a  1009E 
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Fig.  2(g).  Case  A:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  80%  -  RH  <  100% 

RH  2  100% 
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Fig.  2(h).  Case  A:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  so%  s  rh  <  100%  |5gi#i| 

nn  2  100% 
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Fig.  2(1).  Case  A:  analysis  of  observed  surface  Isotherms  (*K)  and  derived 
50-m  wind  field  at  verification  time  (arrows  Indicate  wind  direction;  numbers 
Indicate  wind  speed  In  m  sec-1). 
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Fig.  2(j).  Case  A:  prognosis  of  surface  isotherms  (°K)  and  50-m  wind  field 
valid  at  verification  time  (arrows  indicate  wind  direction;  numbers  indicate 
wind  speed  in  m  sec-*). 
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Fig.  2(k).  Case  A:  vertical  cross  section  along  grid  row  M  =  4  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time,  ho','  rh  <  ioo%Pjgl  rh  *  100%  KSU 


Fig.  2(1).  Case  A:  vertical  cross  section  along  grid  row  M  =  5  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  hW  mi  <  ioo%i|gg|  rh  2  1003, 
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Fig.  2(n).  Case  A:  vertical  cross  section  along  grid  row  M  =  7  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  nor*'  s  rh  <  100%  fill  rh  a  100% 


Fig.  3.  Case  A:  analysis  of  observed  surface  isotherms  (°K)  and  derived 
50-m  wind  field  at  initial  time  (arrows  indicate  wind  direction;  numbers  indicate 
wind  speed  in  m  sec"1). 
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Fig.  4.  Case  A:  predicted  terrain-induced  vertical  velocity  (cm  sec 
verification  time. 


)  at 
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Fig.  5.  Case  A:  net  sensible  heat  (Cal  cm 
during  the  forecast  period. 


)  added  to  the  transition  layer 
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Fig.  6.  Case  A:  net  latent  heat  (Cal  cm 
during  the  forecast  period. 


)  added  to  the  transition  layer 


— 
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Fig.  7(a).  Case  B,  observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  initial  time.  Shading  indicates  regions  covered  by  broken  or 
overcast  low  cloudiness  ,  or  experiencing  precipitation  kill  . 
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Fig.  7(b).  Case  B,  observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  verification  time.  Shading  indicates  regions  covered  by  broken 
or  overcast  low  cloudiness  ,  or  experiencing  precipitation  FsM  . 
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Fig,  7(c).  Case  B:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80%  s  rh  <  ioo%  j&pgi 
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Fig.  7(d).  Case  B:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  80%  *  rh  <  100% 


54 


Fig.  7(e).  Case  B:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  80%  ^  RH  <  ioo% 
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Fig.  7(f).  Case  B:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80%  ^  rh  <  100% 

rh  *  100% 
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Fig-  7(g).  Case  B:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  »a%  s  rh  <  100% 

mi  2  ioo» 
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Fig.  7(h).  Case  B:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  80%  s  rh  <  100%  [fill 

RII  2  1005 
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Fig.  7(i).  Case  B:  analysis  of  observed  surface  isotherms  (°K)  and  derived 
50-m  wind  field  at  verification  time  (arrows  indicate  wind  direction;  numbers 
indicate  wind  speed  in  m  sec"1). 
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Fig.  7 (j) .  Case  B:  prognosis  of  surface  isotherms  (°K)  and  50-m  wind  field 
valid  at  verification  time  (arrows  indicate  wind  direction;  numbers  indicate 
wind  speed  in  m  sec"*). 
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Fig.  7(k).  Case  B:  vertical  cross  section  alon^  grid  row  M  =  -1  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  sn<  mi  <  100%  iff}  Wi  -  100%  ^3 
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Fig.  7(1).  Case  B:  vertical  cross  section  along  grid  row  M  =  5  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  sW  *-*  iui  <  100%  f||||  rh  2  100% 


Fig.  7(m).  Case  B:  vertical  cross  section  along  grid  row  M  =  6  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  ^  rh  <  100%  ESggl  rh  2  100?, 
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Fig.  7(n).  Case  B:  vertical  cross  section  along  grid  row  M  =  7  showing  predicted  temperature  (°K)  and  humidity, 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  ro%  srh<  MKftHU  rh  a  100%  R53 
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Fig.  8.  Case  B:  south-to-north  vertical  cross  section  of  temperature  (°K)  and  humidity  along  the  coordinate  line 
5  (solid  line  indicates  frontal  position),  no^f-  £  RH  <  ioo%  fej  rh  a  ioo%  RSM 


Height,  meters 


Fig.  9.  Case  B:  initial  and  predicted  temperature  structure  at  grid  point 
L  =  9,  M  =  6,  and  observed  (at  verification  time)  structure  at  a  nearby  radio¬ 
sonde  station. 
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Height,  meters 


Fig.  10.  Case  B:  initial  and  predicted  temperature  structure  at  grid  point 
L  =  6,  M  =  6,  and  observed  (at  verification  time)  structure  at  a  nearby  radio¬ 
sonde  station. 
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Height,  meters 


Fig.  11.  Case  B:  initial  and  predicted  temperature  structure  at  grid  point 
L  =  5,  M  =  7,  and  observed  (at  verification  time)  structure  at  a  nearby  radio¬ 
sonde  station. 
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Height,  meters 


-10  0  10  20  30 


Temperature,  °C 

Fig.  12.  Case  B:  initial  and  predicted  temperature  structure  at  grid  point 
L  =  7,  M  =  7,  and  observed  (at  verification  time)  structure  at  a  nearby  radio¬ 
sonde  station. 
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Height,  meters 


Fig.  13.  Case  B:  initial  and  predicted  temperature  structure  at  grid  point 
L  =  2,  M  =  8,  and  observed  (at  verification  time)  structure  at  a  nearby  radio¬ 
sonde  station. 


70 


Height,  meters 


Fig.  14.  Case  B:  temperature  and  dew-point  observations  at  radiosonde 
station  226  for  initial  (fine  lines)  and  final  (heavy  lines)  times  of  the  forecast,  and 
their  values  (medium  lines)  predicted  at  the  nearby  grid  point. 
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Fig.  15(a).  Case  C,  observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  initial  time.  Shading  indicates  regions  covered  by  broken  or 
overcast  low  cloudiness  ,  or  experiencing  precipitation  . 
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Fig.  15(b).  Case  C,  observed  data:  surface  isobars  (mb)  and  700-mb  contours 
(tens  of  meters)  at  verification  time.  Shading  indicates  regions  covered  by  broken 
or  overcast  low  cloudiness  ,  or  experiencing  precipitation  . 
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Fig.  15(c).  Case  C:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80%  £  rh  <  ioorr 
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Fig.  15(d).  CaBe  C:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  80%  tut  <  io</V  B 
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Fig.  15(e).  Case  C:  analysis  of  observed  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  80%  ^  kh  <  ioorr  ($:|| 
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Fig.  15(f).  Case  C:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  500  m  above  terrain  height.  80%  s  rh  <  ioo1;  tggga 

mi  a  io<« 
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Fig.  15(g).  Case  C:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  1150  m  above  terrain  height.  80%  s  rh  <  100%  pm 

nn  2  100% 
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Fig.  15(h).  Case  C:  analysis  of  predicted  temperature  (°K)  and  humidity  at 
verification  time  for  level  2000  m  above  terrain  height.  go%  s  rh  <  ioor, 

mi  ^  too'* 
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Fig.  15(i).  Case  C:  analysis  of  observed  surface  isotherms  (°K)  and  derived 
50-m  wind  field  at  verification  time  (arrows  indicate  wind  direction;  numbers 
indicate  wind  speed  in  m  sec"1). 
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Fig.  15(j).  Case  C:  prognosis  of  surface  isotherms  (eK)  and  50-m  wind  field 
valid  at  verification  time  (arrows  indicate  wind  direction;  numbers  indicate 
wind  speed  in  m  sec-1). 
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Fig.  15(1).  Case  C:  vertical  cross  section  along  grid  row  M  -  5  showing  predicted  temperature  (°K)  and  humidity 
and  observed  sky  condition  and  weather  (at  top  of  figure),  valid  at  verification  time.  h<v;  nn  <  rh  2  100% 


r3^ 


•2  0) 


Height,  meters 


Fig.  16.  Case  C:  observed  initial-time  and  predicted  temperature  at  grid- 
point  L  =  10,  M  =  7,  and  observed  verification-time  temperature  at  radiosonde 
station  518. 
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regions  reporting  broken  or  overcast  low  cloudiness  ( i .  i 
at  verification  time.  Type  of  precipitation  is  indicated. 


Fig.  17.  Case  C: 
and  precipitation^! 


Fig.  18.  Case  C: 
2000-m  levels. 


composite  of  predicted  humidity  for  the  500- , 


1150-,  and 
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Fig.  19.  Case  C:  composite  of  observed  humidity  for  the  500- ,  1150-,  and 
2000-m  levels  at  verification  time. 
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Height,  meters 


Fig.  20.  Case  C:  initial  and  predicted  temperature  at  grid  point  L  =  6, 

M  =  6,  and  observed  verification-time  temperature  at  radiosonde  station  425. 
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Fig.  21.  Case  C:  initial  and  predicted  temperature  at  grid  point  L  =  5, 

M  =  3,  and  observed  verification-time  temperature  at  radiosonde  station  311. 
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Height,  meters 


Fig.  22.  Case  C:  initial  and  predicted  temperature  at  grid  point  L  =  5, 

M  =  7,  and  observed  verification-time  temperature  at  radiosonde  station  429. 
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Height,  meters 


Fig.  23.  Case  C:  initial  and  predicted  temperature  at  grid  point  L  =  7, 

M  =  7,  and  observed  verification-time  temperature  at  radiosonde  station  520. 
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Fig.  24.  Case  C:  initial  and  predicted  temperature  at  grid  point  L  =  7, 

M  =  4,  and  observed  verification-time  temperature  at  radiosonde  station  317. 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 


Based  on  the  limited  number  of  tests  conducted,  we  may  tentatively  conclude  that 
the  model  is  capable  of  producing  realistic  and  fairly  accurate  predictions  of  the  structure 
of  the  atmospheric  boundary  layer.  Our  recommendations  for  further  research  using 
this  or  similar  models  are: 

(a)  Attempt  to  improve  the  strong-inversion  formulation  for  the  contact 
layer  equations,  especially  in  situations  with  warm  advection. 

(b)  Formulate  a  procedure  for  incorporating  the  computation  of  the 
evaporation  of  large-scale  precipitation  within  the  boundary  layer,  and  for 

the  computation  of  precipitation  from  super-saturated  portions  of  the  boundary 
layer. 

(c)  Investigate  the  possibility  of  computing  infrared  radiative  flux 
divergence  for  cloudy  situations  within  the  boundary  layer. 

(d)  Examine  the  potential  for  developing  improved  diagnostic  formulas 
for  the  horizontal  wind,  including  the  dependence  of  the  eddy  viscosity  on  the 
stability  of  the  air. 

(e)  Develop  objective  analysis  techniques  for  combining  surface  and 
radiosonde  observations  into  grid-point  data  within  the  boundary  layer. 

It  is,  of  course,  desirable  that  a  continuous  effort  be  made  to  interpret  the  results 
of  fundamental  theoretical  and  empirical  investigations  of  the  structure  of  the  boundary 
layer  in  the  light  of  the  requirements  for  synoptic-scale  weather  forecasts. 
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APPENDIX 


LOGICAL  AND  COMPUTATIONAL  FORMULATION  OF  THE 
NUMERICAL  MODEL 

The  equations  of  the  boundary  layer  model  were  converted  to  a  form  suitable  for 
solution  on  a  digital  computer.  A  set  of  specifications  for  a  computer  program  was 
prepared  to  guide  the  programmer  in  writing  the  code.  The  actual  program  is  avail¬ 
able  at  the  United  Aircraft  Research  Laboratory,  East  Hartford,  Connecticut  or  through 
request  to  the  433L  Systems  Program  Office,  Electronics  Systems  Division,  Air  Force 
Systems  Command,  L.  G.  Hanscom  Field,  Bedford,  Mass.  Including  extensive  output 
at  hourly  intervals,  the  time  required  for  a  12-hr  forecast  on  a  7094  computer  was 
27  minutes. 

In  this  Appendix,  we  present  the  logical  and  computational  basis  for  the  numerical 
model. 

Data 

Input  data  will  be  provided  in  tabulated  form.  The  format  for  the  machine  input 
is  left  to  the  discretion  of  the  programmer. 

In  the  table  below,  we  indicate  the  symbol,  name,  system  of  units,  approximate 
significances  and  magnitude,  and  storage  requirements  for  the  several  variables  used  in 
the  computation.  Those  quantities  to  be  provided  as  input  are  marked  with  an  asterisk. 
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TABLE  VI 
DEFINITIONS 


Symbol 

Definition 

Units 

Magnitude  and 
significance 

Storage 

T 

Air  temperature* 

•c 

XXX  .X 

10-10-12-2 

q 

Specific  humidity* 

none 

XX.X-  Ex 

10-10-12-2 

r 

Specific  moisture* 

none 

XXX*  Ex 

10*10*12-2 

UGH1 

x— component  of  geostrophic  wind 
at  Z=H  at  Initial  time* 

cm  sec'"1 

XX.X-  Ex 

10-10 

UGH2 

x~componcnt  of  geostrophic  wind 
at  Z-H  at  +6  hrs* 

cm  sec-1 

XXX-  Ex 

10-10 

UGH3 

x-component  of  geostrophic  w  ind 
at  Z=H  at  +12  hrs* 

cm  sec”1 

XX  .X-  Ex 

10-10 

VGH1 

y-component  of  geostrophic  wind 
at  Z=H  at  an  initial  time* 

cm  sec"1 

XX.X*  Ex 

10*10 

VGH2 

y-component  of  geostrophic  wind 
at  +6  hrs* 

cm  sec"1 

XXX-  Ex 

10*10 

VGH3 

y— component  of  geostrophic  wind 
at  Z^H  at  +12  hrs* 

cm  sec"1 

XXX-  Ex 

10*10 

ICLU1 

Upper  cloud  index  0— +6  hr* 

none 

X 

10*10 

ICLU2 

Upper  cloud  index  +6— +12  hr* 

none 

X 

10-10 

WET 

Ground  moisture  factor* 

none 

XXX 

10-10 

E 

Height  of  terrain* 

cm 

XXXX-  Ex 

10-10 

zo 

Surface  roughness  factor* 

cm 

XXX*  Ex 

10-10 

XM 

Map  scale  factor* 

none 

X.XXXX 

10-10 

TS 

Surface  temperature* 

•A 

XXXX 

10-10 

RHS 

Surface  relative  humidity* 

none 

XXX 

10-10 

f 

Coriolis  parameter* 

sec"1 

XXXX-  Ex 

10*10 

P 

Mean  value  air  density* 

gm  cm"^ 

XXXX-  Ex 

1 

(°Z)k 

Vertical  grid  spacing* 

cm 

XXX .XX*  Ex 

13 

Cl 

Sine  of  solar  declination 

none 

X.XXXXXX 

1 

C2 

Cosine  of  solar  declination* 

none 

X.XXXXXX 

1 

C3 

Tangent  of  solar  declination* 

none 

X.XXXXXX 

1 

D 

Standard  horizontal  grid  interval* 

cm 

XXX. XX-  Ex 

1 

DT 

Time  step  interval* 

sec 

XXXXX 

1 

TIME  0 

Initial  time  (EST)* 

sec 

XXXXX.X 

1 

K 

Gravity  acceleration* 

-9 

cm  sec  “ 

XXXX 

1 

C 

P 

Specific  heat* 

ergs  gm"1  deg"1 

XXXXXXX-Ex 

1 

y 

Ratio  of  g  to  Cp* 

deg  cm"1 

XXXX-Ex 

1 

0 

Mean  air  temperature* 

-A 

XXXX 

1 

& 

Forced  convection  parameter* 

none 

XXX 

\ 

Free  convection  parameter* 

none 

X.XX 

1 

K, 

Von  Kar man's  constant* 

none 

X.XX 

1 

zi 

Instrument  shelter  height* 

cm 

XXX. 

10*10 

SLAT 

Latitude* 

radians 

XXXXXXX 

10-10 

SLNG 

Longitude* 

degrees 

XXXXX 

10-10 

♦To  be  provided  as  input. 
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TABLE  VI  (Continued) 


Symbol 

Definition 

Units 

Magnitude  and 
significance 

Storage 

ItA 

Radiation  coefficient* 

°F  hr"1 

.XXXX  Ex 

10-10-3 

RB 

Radiation  coefficient* 

“F  hr-1 

.XXXX  •  Ex 

10*10-3 

RC 

Radiation  coefficient* 

°F  hr*1 

.XXXX  •  Ex 

10-10-3 

RD 

Radiation  coefficient* 

*F  hr"1 

.XXXX  •  Ex 

10-10-3 

K 

Mixing  coefficient 

cm2  sec-1 

10-10-12 

K 

SFC  mixing  coefficient 

cm2  sec"1 

10*10 

w 

Frictional  vertical  velocity 

cm  sec"1 

101013 

w 

Terrain  vertical  velocity 

cm  sec'1 

10-10*13 

u 

x-component  horizontal  wind 

cm  sec'1 

10-10-13 

V 

y-component  horizontal  wind 

cm  sec'1 

10-10-13 

i 

Convection  index 

none 

10-10 

RATIO 

Wind  speed  ratio 

none 

10-10 

VF 

Vapor  flux  (U,Q*) 

cm  sec-1 

10-10 

HF 

Heat  flux  (U*0J 

deg  cm  sec'1 

10-10 

Friction  velocity 

cm  sec'1 

10-10 

<P 

Geostrophic  deviation  angle 

none 

10-10 

q 

SFC  specific  humidity 

none 

10-10 

zk 

Height  of  k1*1  level 

cm 

13 

a 

Gaussian  elimination  coefficients 

none 

13 

b 

Gaussian  elimination  coefficients 

none 

13 

c 

Gaussian  elimination  coefficients 

none 

13 

TIME  0 

Initial  time 

sec 

1 

TIME 

Time 

sec 

1 

TIMER 

Radiation  time 

hr 

1 

S 

Wind  speed  at  Z  =  h 

cm  sec'1 

10-10 

U 

x-component  wind  at  Z  =  h 

cm  sec-1 

10-10 

V 

y-component  wind  at  Z  =  h 

cm  sec"1 

10-10 

ui' 

x-component  geostrophic  wind 
at  Z  =  H 

cm  sec"1 

10-10 

v£ 

y-component  geostrophic  wind 
at  Z  =  II 

cm  sec-1 

10-10 

ug 

x-component  geostrophic  wind 

cm  sec'1 

Share  with 

Vg 

y-component  geostrophic  wind 

era  sec'1 

U  and  V 

A 

Dummy  field 

10-10 

B 

Dummy  field 

10-10 

C 

Dummy  fleid 

10*10 

R 

Dummy  field 

10-10 

AA 

Dummy  field 

10-10-13 

DUMYA 

Dummy  field 

10-10 

DUMYB 

Dummy  field 

10-10 

♦To  be  provided  as  input. 
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Computations 

Coordinate  System  and  Finite  Difference  Notation 

The  independent  variables  are  the  three  space  dimensions  (x,y,  z)  and  time  (t). 

The  x,  y  coordinates  are  the  horizontal  coordinates  (approximately  parallel  to  the  earth’s 
surface).  The  z  coordinate  is  directed  normal  to  the  earth’s  surface  and  counted  positive 
toward  the  zenith. 

We  define  a  finite  set  of  points  in  space  as  the  finite-difference  grid.  The  points 

forming  this  grid  will  have  the  coordinates,  X^,  Y  ,  Z^,  with  £  =  1,  2,  ...,  L;  m  =  1,  2, 

...,  M;  k  =  1,  2,  ...,  K.  Thus,  there  will  be  L  M  K  grid  points  covering  the  space  within 

which  our  equations  are  to  be  solved.  We  will  compute  the  solution  only  at  discrete  times, 

t  ;  n  =  1,  2,  ....  N. 
n 

The  intervals  between  grid  points  are  defined  as  follows: 


<DX»t  - 

Xl' 

xn 

£  =  2,  ...,  L 

<DY>m  “ 

Y 

m 

1 

*< 

5 

1 

h- 1 

m  =  2,  ...,  M 

<DZ>k  ' 

V 

zk-i 

k  =  2,  ....  K 

<DT,n  = 

t  - 
n 

ln-l 

n  =  2,  ...,  N 

In  this  work,  we  will  take, 

D  =  (DX)^  =  (DY)^  (constant) 

DT  =  (DT)^  (constant) 

If  a  dependent  variable  is  defined  [e.g.,  F  =  F(x,y,z,t)],  we  will  use  the  following  notation: 


F*'"  ‘  F 

£ ,  m 


F(x  =  x  ,  y  =  y  ,  Z  =  Z  ,  t  =  t  )  . 
£  m  k  n 


If  the  variable  depends  upon  fewer  than  four  independent  variables,  only  the  indices 

required  will  be  used  [e.g.,  F^  =  F(x  =  x.,  y  =  y  ,  t  =  t  )]. 

*  y  m  i  m  n 

Definition  of  an  Advection  Operator 

In  setting  up  the  difference  form  of  the  differential  equations,  we  must  treat  advection 

terms  of  the  form 

3T  ,  3T 

u  ttt  and  v  ~ 
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The  method  of  upwind  differencing  will  be  applied  to  these  terms  (see  [9]).  We  will  use 
k  k 

the  symbols  [X(T)]^  m  and  (Y(T)]^  ^  to  represent  these  terms.  The  symbols  are 
defined  as  follows: 

For  f  =  1  and  m  =  1 : 

k  k 

If  u  >  0,  or  v  >  0,  set 

i. ,  J.  X  ,  X 

[X(T)]k  ^  =  [Y(T)I^*”  =  0 
Otherwise,  set 


(X(T)J^  J  = 

k 

u  * 

1,1 

XM1  1 
- ±JL-Jl  * 

D 

,„,k,  n  k,  n, 

IT  -  T  1 

2,1  1,1J 

'Y<T>(,  j  = 

k  * 

v  * 

1.1 

XM1  1 
- ij-i  * 

D 

_k,n  ^k,n, 

|ti,  2  *  Ti:i* 

For  £.  -  L  and  m  =  1 : 

If  u^  <  0,  or 

L,  1 

k 

v  >0,  set 

L,  1 

Ix<T»kLil  = 

|Y,T)1L,1  ° 

0 

Otherwise,  set 

k 

x_^  , 

D 

(Tk'n  -  Tk’"  ] 
L,  1  L-l,  1* 

1  ■ 

u  * 

L,  1 

'Y<T»L,  1  = 

Vk 

L,  1 

^  . 

D 

k,n  _  k,  n 

L,  2  L,  1J 

For  S.  =  1  and  m  =  M 

• 

If  \  M  '  °>  °r 

Vl,  M  <  °’ 

set 

'X<T>lk  M  " 

0. 

Otherwise,  set 

- 

k  * 

u  * 

1,  M 

XM,  «« 

_ LJI  * 

D 

k,  n  k,  n 

1  2,M  1  1,  M 

lY<T»l  M  - 

k 

v  * 

1,  M 

XM  A/I 
_ L-M  * 

D 

k,n  k,n 

1  1,  M  1,M- 
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For  £  =  L  and  m  =  M : 


If  UL,M  <  °'  0r  VL,M  •=  °'  Set 

Hm  ‘  |y<t<,m  ■  »• 

Otherwise,  set 


[X(T)iL,M 

k 

UL,M 

XM 

*  L.M  + 

D 

(Tk'n 
L,  M 

-  Tk,n  1 
L-1,MJ 

|Y(T)1l,m 

k 

=  V 

L,  M 

XM  .. 

*  ki.M  + 

D 

r_k,  11 
[T 

1  L,  M 

-  Tk’n  ] 

L,  M-1J 

For  £  =  2,  L  -  1 

and  m  = 

2,  ....  M  -  1: 

If  u  —  0 ,  set 

£,  m 


XM 


k,  n 

l ,  m 


D 


fTk,  n  k,  n 

1  i,m  £-1, m  ’ 


If  u.  <0,  set 
£,  m 


[X(T)]k  =  u" 

1  v  /J£,  m  £,  m 


XM 


k,  n 
£,  m 


D 


*  rrpk>  ^  T1^’  ^  1 

1  £+l,m  "  l£,  m  ’ 


If  v 


£,  m 


0,  set 


k  k 

[Y(T)L  =  v* 

£,  m  £,  m 


XM 


k,  n 

l.  m 


D 


k,n  k,n 
1  l.m-11' 


If  V  <  0,  set 
x,  m 


[Y(T)]k  =  vk 

£,  m  £,  m 

For  £  =  1;  m  =  2,  M  -  1; 

If  u  >  0,  set 

1,  m 


XM 


k,  n 

£,  m 


[T 


k,n 

£,  m+1 


Tk’nl 
1  £.,  m  • 


[X(T)]k  m  =  [Y(T)]k  =  0 

l,m  1,  m 


104 


Otherwise,  set 


k 

XT  I  =  u 

1,  m  1,  m 


XM 


1,  m 


D 


IrT,k,n  „,k,n  . 

[T  ’  ‘  Ti\J. 

*_  j  ni  X  j  m 


and,  if  v,  ^  0,  set 
1,  m 


[Y(T)]k  m  =  v!< 

1 ,  m  1 ,  m 


XM. 


1,  m 


D 


rrT,k,  n  mk,  n  , 

IT.’  -  T  ’  ], 

1,  m  1,  m-l 


If  v,  <  0,  set 

1,  m 

k  k 

[Y(T)]‘  =  v, 

l,m  l,m 


XM 


1,  m 


D 


k,  n  Tk>  n  I 

1  1,  m+1  1,  mJ  ’ 


For  £  =  L;  m  =  2,  ...,  M-l: 


If  u„  <0,  set 
L,  m 


|Xm)kL,m  -  [Y(T)]kUm  =  0 


Otherwise,  set 

[X(T))k  =  u* 

m  i-i,  m 

,  k 

and,  if  v  0,  set 

L,  m 

iv(T)]1;  =  k 

JLi,  m  m 


XM 


L,  m 


D 


XM 


L,  m 


D 


,  k,  n 

[T 

L,  m 


k,  n  , 
T 

L-l,  mJ 


IT 


k,n  _  Tk,n 
L,  m  L,  m-l 


If  vT  <  0,  set 
L,  m 


IY(T)]J  =  vk 

L,  m  L,  m 


XM 


L.  m 


D 


rTk»n  -  Tk’n  1 
1  L,  m+1  L,m 


For  m  =  1,  £  =  2,  ...,  L-l: 


If  v£  1  >  0,  set 


lX(T)]£,l  =  lY(T)]£,l  =  0 

Otherwise,  set 


l  Y(T)  Jf  i  =  V£  1 


XM, 


D 


_k,n  ^k,n, 
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and,  if  uk  ^  s  0,  set 


k  1c 

[X(T)L  =  u 

ni,i  t,i 


k 


If  1  <  0,  set 


XM 


1x1 


D 


|Tk'n-  Tk'n  , 

1  1,1  1-1,  1J’ 


[X(T)f  =  uk  * 

,J1,1  1,1 

For  m  =  M;  £  =  2,  ...,  L  -  1  ; 

If  V£,  M  <  °*  Set 

W»f,M  ■  >Y<T»fkM  ’  0 

Otherwise,  set 


XM 


11 


D 


[Tk’n  -  Tk’ni 
1  £+1,1  £,1J- 


[Y(T)]k  =  vk 

k  ;J£,M  £,M 


XM 


£,M 


D 


[Tk,n  -  Tk,n  1 
1  £,M  £,M-1J’ 


and,  if  u.  >  0,  set 
£,  M 

|XKm  ■  <M 

If  ukM  <  0,  set 


XM 


£.M 


[Tk>  n  -  Tk’  n  ] 
£,M  1-1,  M  ’ 


[X(T)J 


k 


=  u. 


£,  M  £,  M 


XM 


£,  M 


D 


[Tk’n  -  Tk>n] 
1  1+1,  M  £,MJ 
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Fig.  25.  Major  computational  steps  for  numerical  model. 
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Fig.  25.  (continued) 
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Fig.  25.  (continued) 
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Read  Input  Data 

The  information  relating  to  input  data  iB  given  above.  We  will  indicate  here 
that  there  are  two  types  of  input  quantities.  The  first  type  consists  of  a  simple  constant 
which  will  be  unmodified  during  the  computation.  The  second  type  consists  of  initial 
values  of  the  dependent  variables  that  are  modified  by  the  computation. 

A  computation  is  to  be  made  which  involves  derivation  of  the  height  of  the  various 
levels  above  the  ground. 

\  =  (DZ)X 


For  k  =  2,  ....  K 


z,  =  z,  . 
k  k-1 


(DZ), 


The  input  values  of  temperature  are  in  degrees  centigrade  and  they  are  converted 

to  degrees  absolute.  The  coefficients,  RA**  ,  which  are  given  as  input  must  be  adjusted 

£,  m 

as  follows: 


RA 


II 

£,  m 


RA 


II 

£,m 


RB**  {1.8 
£,  m  l 


(TS.  -  273.2)  +  32.0} 
I,  m  J 


Printout  Input  Data 

This  output  is  required  to  insure  that  the  appropriate  numerical  values  have  been 
read  into  the  computer  storage.  This  output  must  include  all  the  quantities  indicated  in 
the  input  data  table  by  asterisks.  It  should  also  include  Holerith  characters  sufficient  to 
identify  the  case  for  which  the  data  apply. 

Geostrophic  Wind 

The  computational  formulas  outlined  here  are  used  to  determine  the  values  of  the 
geostrophic  wind  components  which  will  subsequently  be  used  to  derive  the  actual  wind 
components.  These  two  wind  fields  may  share  memory  locations.  For  clarity,  we  will 
use  a  subscript  g  on  the  symbols  for  the  geostrophic  components  in  these  formulas  but 
these  may  be  omitted  in  the  program. 

The  first  step  in  the  computation  involves  the  evaluation  of  the  geostrophic  wind 
shear. 

For  £  =  1,  ....  L;  and  m  =  1,  ...,  M: 
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For  £  = 


For  £  = 

For  m 

For  m  - 

For  £  = 

Finally, 
For  £  = 


0.5 


Z,  -  Z. 
k  1 


k=K 

E 


k=2 


1.0 

„,k,  n 
T. 

_  £,  m 


1.0 

_k-l,n 
T.  J 
I,  m 


2,  L  -  1;  m  =  2,  M  -  1: 

g*XM.  m 

T»  _  - Zi-Hi.  .  *  (A  _  A  1 

£,  m  2.0*D*f  1  £,m+l  £,m-lJ 

Jt,  m 

g*XM  m 

r  =  -  — - L.rr.1 —  *  fA  -  a  i 

£,m  2.0*  D*  f.  1  £+l,m  £-l,mJ 

Jt,  m 


B. 

B 

£,1 

4,2 

B„ 

B„ 

£,M 

£,  M-l 

c. 

=  Cn 

1,  m 

2,  m 

C, 

=  cT 

L,  m 

L-l,  m 

....  M 

-  1: 

B, 

=  Bn 

1,  m 

2,m 

Bt 

=  bt  , 

L,  m 

L-l,  m 

• « * ,  L 

1: 

C„  = 

C. 

£.1 

£,  2 

C„ 

:  c. 

£,M 

£,  M-l 

we  may  compute  the  geostrophic  wind  field. 


L;  m  =  1,  ....  M;  k  =  1, 


K: 
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U6«,m  =  <Ug\?m  *U.0  + 


,l.n 


rTl.n  .  TK.nn 

£,m  £,  m 

_K,  n 
T 

£,  m 


ZK  ~  Zk 

ZK-Z1 


-  t:,_  *  b,  *  |z„  -  z. ) 

£,  m  X,  m  K  kJ 


z  - 

Z,  \ 

K 

k 

z„  - 

Z, 

K 

1/ 

Tl,n  _  JK.n. 

£.m  l.  ml  g. 


£,  m 


rK,n 

£,  m 


IE 


-  E, 


f  2.0  *  D  1  £,m+l  £,  m-1 


,1 


,  n  .  H  n  .  J  ,  A 
=  (v  )  „  *  ^1.0  + 


v  7 "  =  (v")  „ 
g  £,m  g'  £,m 


-Tl,n  _  TK,n-, 
£,  m  i,  m 


£,  m 


ZK  -  Zk 


z  -  z 

K  1 


,l,n 


-  Tt',m  *  °«,m  *  |ZK  -  Zk‘ 


Z  -  Z  \ 
K  k 

T1* n  _  xK,n\ 
£,  m  £.m 

|ZK  ~  Zlj 

TK’n 

'  Tl,m  ' 

Neutral  Stress  and  Geostrophic  Deviation 

(XM)„ 

*  Ja  *  - m  *  ITT  _  £  1 

f  2.0  *  D  lt£+l,m  £-1,  m 


The  stress  is  symbolized  by  u*  and  requires  for  its  computation  the  value  of  the 
geostrophic  wind  speed  at  the  surface.  The  geostrophic  deviation  is  the  angle  (in  radian 
measure)  between  the  surface  geostrophic  wind  and  the  actual  wind  in  the  contact  layer. 

For  £  =  1,  ...,  L  and  m  =  1,  ...,  M: 


/  H\  n  *  J 

=  (U  )  „  *  ^1 , 


A£,  m  ^Ug  ^  £,  m 


0  + 


■  l,n 
T  -  T 
£.  m  £.  m 

K,  n 


K 


Z  -  Z 
K  1 


£,  m 
K.n, 


Z  „  "1 

3|C 

K 

z„  -  z. 

J 

K  1 

1 


rml,n  _ 

rT  -  T 
£,  m  £,  m 

^K.n 

T. 

£ ,  m 


£ 


(XM) 


£,  m 


IE 


-  E, 


f£  m  2.0  *D  £,  m+1  £,  m-1 


J 


-  T1,n  *  B  *  Z 
£,  m  £,  m  K 
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r 


B 


£,  m 


(vH)  n  *|l.0  + 
g  f,m  \ 


"T1,n  -  TK,n 
£,  m  £,  m 


rK,n 

£,m 


K 


-n 


ZK-  Z1 


'T'^’  n  *  P  *7 

£,  m  £,  m  K 


z , 

"T.1,n  -  Tf’n_ 

K 

£,  m  £,  m 

z„  -  z 

K,  n 

K  1 

T 

£,  m  J 

:ff  m  =  I(A.  )  ** 

£,  m  £,  m 

2  +  (B  )  ** 

v  £,  nr 

=  A,  -r  C 

£,  m  £,  m 

£,  n 

=  B,  -r  C 

£,  m  £,  m 

£,  n 

(XM)( 

rr  '  'C  m 

-£ - ^  [E, 


-  E, 


[Note  that  if  C.  is  zero  replace  it  by  100.0] 

£,  m 

R.  =  C,  -r  [f„  *  Z„  .  ] 

£,  m  £,  m  £,  m  0  £,  m 

Rf  m  =  [R  ]  -r  £n  (10.0) 

jt,  m  1,  m 

u*  =  C„  *  [0.07625  -  0.00625  *  R„  ] 
*  £,  m  £,  m 

£,  m 


f£  m  2.0  *D  1  £+1,  m  £-1,  m 

0,1/2 


<Pi,  m 

C„ 

£,  m 

R„ 


3.14159 

180.0 

*  10.625  * 

(R.  )  ** 

£ ,  m 

2  - 

12.75 

=  A„  * 

cos  '’W 

-  Bo  m  * 

sin 

[  <Pt  m 
£,  m 

=  B.  * 

£,  m 

°os  Wt  m) 

+  A„  * 

£,  m 

sin 

[<Pt  m. 
ii,  m 

£,  m 


[Note  that  <p  is  in  radians] 


Contact-layer  Static  Stability 

In  order  to  determine  the  appropriate  value  for  the  stress  and  to  decide  if  free-  or 
forced-convection  formulas  are  to  be  used  in  the  contact  layer,  we  must  compute  the 
static  stability  of  the  contact  layer.  We  therefore  evaluate 

For  £  =  1,  ....  L  and  m  =  1,  ...,  M: 


\m  ‘  lT«,m  -  TC]  +  »  ’  IZ!  -  <Vm  +  ZU,m» 
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Non -neutral  Stress 


Based  on  the  sign  of  A  computed  in  the  equation  above,  we  use  one  or  another 
formula  for  the  stress  (u+) 

For  £  =  1,  L  and  m  =  1,  M: 

If  A  s  0,  set 
£,  m 

u+  =  0.80  *  u+ 

£,  m  £,  m 


u*  =  1-20  *  u+ 

£,  m  £,  m 


The  Limiting  Stability 

Having  computed  the  non-neutral  value  of  u*,  we  may  now  determine  the  limiting 
value  of  the  static  stability  for  the  use  of  forced-convection  formulas. 


For  £  =  1,  ...,  L  and  m  =  1,  ...,  M: 


(u*  )  **  2 

£,  m _ 

4.0*0 


(£n  IZ-l  -  (ZQ  +  Z.)])  **  2 
_ _ _ m  _ 

((K)  **2)  *  [Z1  -  (ZQ  +  Z.)] 


Test  for  Free  or  Forced  Convection 
For  £  =  1,  ...,  L  and  m  =  1,  ...,  M: 


If  A 


£,  m 


>  B 


£,m, 


set 


1 


If  A 


£,  m 


£  B 


£,  m 


£,  m’ 
0 


set 


Free-convection  Formulas 


For  £  =  1,  ... 

,  L 

and  m  =  1,  ...,  M: 

If  lc 
£,  m 

*  0, 

do  nothing. 

If  I, 

=  0, 

set 

£,  m 
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For 


HF 


Jg,  m 


r  2/3  1/3 

^  ^ 

A. 
jg,  m 

3.0  0 

r/ii 

1/3  _ 

1  1 

1/3 

- 

ltzJ 

(Z„  +  Z.) 

'  0,  11 

- 

1/2 


K  =  HF  * 

£,m  jg,m  9 


g  1/3  *  2/3  .  „  4/3 

&  *  \  *  (Z  ) 


£,  m 


u* 

riOO,  +  Z 

jg,  m 

p  3.0  *  u*  *  1.3  -I 

.  £,m 

jg,  m 

m 

K 

\m  J 

-x2/351/3hp.  - 

1  \  1/3 


100 


u. 


1,  m 

122.  +  Z. 


1/3 


WSFC  = 


jg,  m 


in 


l,  m 


0 


jg,  m 


(RATIO) .  =  WSFC  4-  A, 

v  i,  m  £,  m 

B„  =  A,  *  C 

Jg,  m  jE,  m  Jg,  m 

A.  =  A.  *  R, 

jg,  m  £,m  jg,  m 

HF.  =  -  [(HF).  **  3.] 

Jg,  m  u  ;jg,  m  J 

4  6 

Require  10  ^  K.  —  10 

jg.  m 


Forced-convection  Formulas 


-  1, .. 

.,  L  and  m  =  1 

If  I* 
jg,  m 

=  0,  do  nothing. 

If  f/i 

*■  0,  set 

jg,  m 

HF 


jg,  m 


r 9 *in  (Z  -r  (Z  +  Z  ))  *  (u*  )  i 

_ _  i,  m _ i,  m 

2.0  *  g  *  /3  *  *  (Z  -  (Z  +  Z.)) 

JE,  m  1 


4.0  *  g  *  *k  *(Z1  -  (ZQ 


1.0  + 


Jg,  m 


+  Z.))  *  A. 
i  jg,m 


1/2 


0*(u*  *jgn(Z14-(ZQ  +  Z.)))**2 
jg,  m  1,  m 


115 


If  A.  >0,  compute 
£ ,  m 


DUMYA 


-  MM 

1/3  g 


1/2 


u,  2  (Zx  -  Z.  ^  ♦ 

£,  m  £,  m 


£,  m 


1/2 


DUMYB  =  KZ^  4-  (l.O  +  (^g/£Z1HF|  m  4-  (6u,  3  ))] 
£,  m  ’  £,  m 

If  DUMYA  s  DUMYB,  set 

4 

K.  =  10 

£,  m 

I1F  =  10_2*K  *A  v  [Z  -  (Z  +Z  )] 

£,  m  £,  m  £,  m  1  i„  0. 

£,  m  £,  m 

l.n  _ _  l>n  **««  l/2. 


A.  =  0.41  *  (u  ’  *  *  2  +  v 

£,  m  v  g.  g. 

&£,  m  b£,  m 

RATIO  =  2.0  *  ut  4-  A„ 

£,  m  *  £,  m 

’  £,  m 

B  =  A.  *  C . 

£,  m  £,  m  £,  m 

A„  =  A.  *  R, 

£,  m  £,  m  £,  m 


2)  ] 


If  DUMYA  >  DUMYB,  set 

K  =  DUMYB 
£,  m 


£,  m 


u .  £n 

Z,  4-  Zn 

* 

1  o. 

£,  m 

£,  m  J 

Z-Zn  I  HF 
1  0.  I  £,  m 

 £,  m'  


©Uj 


£,  m 


u„ 


£  m 

WSFC  =  £n 


122.  h  Z 


/  Zo 

£,  m  '  £,  m 


RATIO.  =  WSFC  4-  A„ 

£,  m  £,  m 

B.  =  A.  *  C„ 

£,  m  £,  m  £,  m 

A.  =  A.  *  R, 

£,  m  £,  m  £,  m 


If  A,  —  0,  compute 
£,  m 


K  =  kz  ,u* 
£,  m  1  * 


£,  m 


1.0  + 


KZ1  HFlm 

eu! 


116 


For  £ 


If  K.  ^  0,  set  I„  =0  and  use  free  convection. 
£,  m  £,  m 


If  K„  >  0,  compute 

£,  m 


CON  = 


*2/3  IK 


3.0  \0/ 


1/3 


[<zi 


l  I/T 


£,  m 


+  Z; 


J£,  m  l£,  m 


1/3 


1/2 


If  ABSF  (HF£  m)  <  (CON)  ,  set 

I.  =0  and  use  free  convection. 
£,  m 

If  ABSF  (H F  )  ^  (CON)3,  set 

£,  m 


A 


f  m 
1,  m 


Z  /  Z 
l'  0 


£,  m 


£,  m 


P  S 


Z1  '  Z0  |  HF 

x,  m  x,  m 


©Uj 


£,  m 


WSFC  - 


£,  m 


to  ((122.  ♦  Z0  )  *  Z0lim] 


RATIO.  =  WSFC  -r  A 

£,  m  £,  m 

B  =  A  *  n 

£,  m  £,  m  £,  m 

A.  =  A„  *  R 

£,  m  £,  m  £,  m 

*■*  4  6 

[Note  that  in  each  loop  K  is  to  be  forced  to  lie  between  10  and  10  .) 

Horizontal  Wind  in  the  Transition  Layer 

We  recall  here  that  u  (v  )  and  u(v)  will  share  storage  locations. 

S  S 

=  1,  ...,  L  and  m  =  1,  ...,  M: 

C  =  a  -  v1,  ^ 

£,  m  £,  m  g  £,  m 

B.  =  B.  -  u1’  n 

£,  m  £,  m  g  £,  m 


A.  m  =  m  *  <2-0  *  Kc  J] 

£,  m  £,  m  £,  m 


1/2 
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For  £  =  1,  L;  m  =  1,  M;  and  k  =  1,  ....  K: 


k  k 

ix,  =  u  .  + 

£,  m  g  £,  m 


ex»  >-  \m  *  <Zk  -  ZX»] 


)] 


\m  *  C0S  '\m  *  <Zk  ‘  Z1 

Sin  l\m  ’  (Zk  '  Zl»] 


+  C.  * 
£,  m 


V  =  V  + 

£,  m  g  £,  m 


exp  At,m  *  (Zk  '  Zl)l_ 

Km  *  C°S  [AjE, m  ’  <\  '  Zl» 
-  B,  _  *  sin  [A,  _  *  (Zk  -  Z1)j] 


£,  m 


£,  m 


Divergence  of  the  Horizontal  Wind  Field 

This  computation  is  required  as  a  preliminary  to  the  calculation  of  the 
“frictionally-induced”  vertical  velocity. 


For  £  =  2,  L  -  1;  M  =  2,  M  -  1;  k  =  1,  K 

k 


(AA) 


,  XM, 
k  _ l^rn  + 


£,  m  2.0  *D 


^U£+l,m  U£-l,m^  +  ^V£,  m+1  V£,m-1^ 


Frictionally- i ndu ccd  Vertical  Velocity 

The  parameter  is  denoted  by  w  and  is  calculated  as  follows: 
For  £  =  2,  ...,  L  -  1  and  m  =  2,  ...,  M  -  1: 


w,  =0 
£,  m 

w*  =  -  0.5  ♦  (DZ)  *  (AA)  “  +  (AA)  1 


£,  m 


'2 


£,  m 


£,  m 


For  k  =  3,  ....  K: 


w.  k  =  ^  -  0.5  *  (DZ)  * 

Jt,  m  m  k 


(AA)  k  +  (AA)k  1 
v  '£,  m  '£,  m 


If  £  =  1,  L  or  m  =  1,  M,  w,  =  0  for  all  k  . 

f ,  m 
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Terrain-induced  Vertical  Velocity 

The  motion  of  the  air  normal  to  a  terrain  gradient  forces  an  ascent  or  descent  to 
conform  to  the  terrain  surface.  This  vertical  motion  is  computed  here  using  the 
operator  defined  in  Subsection  3.2. 


For  £  =  1,  L;  m  =  1,  ....  M;  k  =  1,  ...,  K: 

w  k  =  [X(E)J  k  +  [Y(E)]  k 
Jt,  m  Jt,  m  id,  m 

Richardson’s  Number  and  the  Mixing  Coefficients  in  the  Transition  Layer 


For  £  =  1,  ....  L;  m  =  1,  ...,  M;  k  =  2,  K: 
„k,  n  _k-l, 

ty  -  t  ■ 

£,  m  i,  m 


A  = 


(DZ), 


+  y 


f  k,n  k-l,n 

P  =  u.  -  u„ 

I  £,m  £,  m 

If  A  >  0, 

and  if  P  =  0,  set 

K  k  =  104 
f,  m 

but  if  P  *  0,  set 


BRI  =  ^  (DZ)2  *  A  -r  P 


k,  n  k-l,n 

v  -  v 
£,  m  £,  m 


If  BRI  >  1.0,  set 


K^'"  =  104 
£,  m 

If  BRI  <  1.0,  set 

(Z, 


K 


k,  n 
£,  m 


Zk  l)  (1  -  BRI) 


2.0 


Pl/2 

(DZk 


If  A  <  0, 

and  if  P  =  0,  set 


K 


k.  n  _  X<Zk  -  2k-l* 
'£,  m  4.0 


SA 

0 


1/2 
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If  1,  m  =  1, 
£,  m 


and  if  K.  =  10  ,  set 
£,  m 


VF.  =  K. 

£,  m  £,  m 


l.n 

q  -  q 

£.  m  £.  m 


’l  ~  ^Z°P  m  +  Zh  rJ 

i,  m  £,  m  J 


but  if  K.  >  10  ,  set 

£,  m 


VF 


£,  m 


-  q 


u*( 


£.  m  £.  m _ £,  m 


£n 


J£,m 


+  Z. 


^«*«*IZi-(Z»  +Zi)l*PF  I 

y £,  m | 

[9  *  (u*  )  **3.] 

£,  m 


Values  of  T,  q,  and  ratZ  =  Hatt  =  t 


n+1 


For  £  =  1 . L  and  m  =  1,  ....  M: 

mK,n+l  _K,  n  * 

T  ’  =  T  ’  +  (DT)  * 

£,  m  £,  m 


-  [X(T)J  K  -  [Y(T)]  K 
£,  m  £,  m 


“  Y 


_  K  K 

W.  +  w« 

£,  m  £,  m 


)  '  *  [ 


K,n  K-l,n 

£.m  £.m 


(DZ) 


K 


K,  n+1  K,  n 


q,’  =  q*  +pt) 

£,  m  jLj  m 


<x<<m  -  rrwitKm 


-  w 


K,  n  _  K-l,n/l 
K.  *  |  ^f,  m  ^£,  m 


r 

K,  n+1 

r. 

£,  m 


£,  m 


=  rK’*  +  (DT)  * 

£,  m 


(DZ) 


-  w 


K 

£,  m 


K 


-  lXMl«Km  - 


K,  n  K-l,  n  q 

r  -  r  ■< 

£,  m  £,  m 


(DZ) 


K 


Gaussian  Elimination  Coefficients  (Temperature,  Specific  Humidity,  and  Specific 

Moisture)^ 

These  four  blocks  are  united  here  since  a  separate  description  of  each  section 
would  be  difficult  and  the  programmer  will  probably  use  the  same  basic  routine  for 
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each  of  the  3  equations. 

The  differential  equations  to  be  solved  are  replaced  by  difference  equations. 

If  we  let  U  represent  a  general  dependent  variable,  the  difference  equations  are 
all  of  the  form 

k  k-l,n+l  ,  k  TTk,n+l  k  ,Tk+l,n+l  ,k 
a.  U.  +  b,  U„  +  c.  U.  =  d, 

£,  m  £,  m  £,  m  £,  m  £,  m  £,  m  £,  m 

(k  =  1 .  K-l) 

This  is  a  special  form  of  the  general  matrix  equation, 


A.  .  U.  =  D.  (i  =  1,  ...,  K-l;  j  =  1,  ...,  K-l) 

1 1 J  J  t 

in  which  A.  .  is  a  (K  -  1)  (K  -  1)  matrix.  In  our  case,  only  the  three  main  diagonals  of 
L  J 


A.  .  are  non-zero. 
i.J 

The  method  for  solving  the  system  is  an  adaptation  of  the  Gaussian  elimination  scheme 
described  in  the  book  by  Forsythe  and  Wasow  [9],  In  the  following  we  will  describe  the 
method. 


Method 

Step  1:  Find  the  element  in  A.  .  which  has  the  greatest  absolute  value. 

J 

Let  the  element  be  A  (i.e.,  the  coefficient  for  which  i  =  R  and  j  =  S). 

lv,  b 


and 


Step  2:  Set 

A  (1)  =  A  —  a 
AR,j  AR,j  ’  AR,S 


=  D  -r  A 

R  R  R,  S 


(j  =  1,  ....  K  -  1) 


and 


Step  3:  For  all  i  (except  i  =  R)  and  for  all  j,  set 

(1) 


A(1)  =  A 
l.J 


-  A.  *  A' 

LJ  i,S  R,  J 


d!1*  .  D.  -  A  *  D*1* 
i  i  i,S  R 


v(l) 


Step  4:  Find  the  element  of  the  modified  matrix  A'  .  which  is  greatest  in 

L  J 

absolute  value;  exclude  from  consideration  all  elements  for  which  i  =  R.  Let  the  element 

(1) 


chosen  be  A 


P,  T 
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Step  5:  Set 

A«  -  A(1)  -  A(1) 

Ap,i  AP,i  •  ap,t 
0<2)  .  D*11  *  *« 


(j  =  1,  K  -  1) 


Step  6:  For  all  i  (except  i  =  P)  and  for  all  j, 


a!2!  =  A*1!  -  a'1'  •  a<2>. 
i,J  l.J  i,T  P,J 

D(2)  =  D(1>  -  A<X>  *  D(2> 
i.j  i.j  i»  T  P 

Now  the  coefficients  a(2).  are  scanned  to  determine  the  element  largest  in  absolute 

i.J 

value  (those  elements  are  excluded,  however,  for  which  i  =  P  or  R), 

The  basic  steps  are  iterated  until  one  arrives  at  the  final  coefficient  matrix 

/K_i\ 

A.  .  .  This  matrix  will  have  just  one  non-zero  element  in  each  row  which  will  have 

i.J 

the  value  unity.  The  solution  for  the  U.  is  found  then  as  follows: 

J 

If  A:  '  is  unity, 

I,  J 

U.-D?"1*. 

J  J 

We  note  that  at  each  horizontal  grid  point,  each  of  the  three  unknowns  must  be 
obtained  by  an  application  of  this  method  using  the  appropriate  coefficients. 


Derivation  of  the  Temperature  Difference  Equations 

If  we  omit  the  term,  (dT/dt)^,  the  differential  equation  may  be  replaced  as  follows: 
For  k  =  1,  £  =  X,  m  =  M : 


T1’n+1  _  t1’0 

4 

Aii _ 1  hiR  _  [xmll 

(DT)  1  T  A.,m 

-  tV(T)]^  - 

-1 
y  w 

\,p 

2,n+l  l,n+l^ 

_ 

2.0 

+  * 

2  r 

K,  *  i  y  + 

LAxiL 

-  hfn 

(DZ)2 

M  |  7 

(DZ)2  j 

X.p 
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Fox'  k  =  2,  ...,  K  -  1;  £  =  A,  m  =  [i: 

Tk,n+1  Tk,n 

-  |DT)  ^  ■  -  |X(T)lXk,-  mT)C  -  +  <M. 

rTk+l,n+l  _  Tk-l,n+l-| 

-X.M  Xju 


-  w 


KDZ>k+l  +  (DZ)k] 


2.0 


[<DZW  +  (DZ)k] 


K*+1  * 


k+1,  n+1  k,n+l 
X.M  ~  \,n 


(DZ) 


k+1 


-  K.  * 


Tk,  n+l  _  Tk-l,n+l 

_Kp- _ 

(DZ). 


k.n+1 


Combining  the  coefficients  of  the  unknowns,  T  ’ 

A,  H 

algebraic  equations  may  be  written 

k  „k-l,n-l  k  k,n+l 


(k  =  1,  K  -  1),  the  system  of 


a.  T.  ’  +  b^  T  '  +  c. 

\,M  X.M  \,n  X,n  \,H  \,M 

in  which  the  coefficients  are  defined  as  follows: 

For  k  =  1 


Tk+l,n+l  _  d  k 


a  =0 


\  1  + 

2.0  *  K.  i 

L  (DT) 

l(DZ)2]2  J 

2.0  *  K  “ 
X.M 

UDZ)2]2 

*  ~  1 
y  *  W  + 

X.M 

vmL1.,  + 

A, 

[y(t)i,;m 

1,  n 

T  ’ 

2 

2.0  *  y  *  K,  2.0 

X, u 

(DT) 

(DZ)2 

(DZ), 
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For  k  =  2,  ....  K 

k 

a. 

\,H 


-  2: 


n  _  k  K\,  n 

0.5  w  +  ~ — u~ 
\,ti  (DZ), 


0.5  *  I(DZ)k+l  +  (DZ)k) 


^  k 

Xm 

i  k 


(DT) 


2.0 


<DZ)k+1  -  (DZ)k 


K 


k+1 

X  M 


K,  ■* 
X.  U 


(DZ)k+l  <DZU 


Kk+1 

-  0.5  w  k  +  - 


\,H  (DZ) 


k+1 J 


0.5  *  [(DZ)k+1  +  (DZ)k] 


T  *  («X  +  \k„|  +  |X(TC  +  |Y(TC 

k,  n 

.  _AJ±  .  _ MIX _  *  I  K  k  -  Kk+1 

(DT)  (DZ)k+1  +  (DZ)k  l  X,M 


For  k  =  K  -  1 : 


K-l 


X,M 


kk_1 

0.5  w^1  +  — ^ji£- 


\,M  (DZ). 


'K-l 


0.5  *  l(DZ)K  +  (DZ)^] 


.  K-l 

b.  =  - 

K-l  „ 

c.  =0 
X,  M 

HK-1 
d  =  y 
X,H 


1 

2.0 

* 

rX  k-1  ,i 

X,U  X.U 

(DT)  + 

«dz)k  +  (dz)k1]J 

[(DZ)k  (DZ)k.JJ 

K-l  ^  K-l 
w.  +  w. 

,K-1,  n 

Xm  ^  _ 


+  + 


2.0  *  y 


(DT)  [(DZ)k  +  (DZ)k_1] 


Kf'1  -  k  K 

\,H  X,M 


+ 

K  1 

L,  K-l 

0.5  w  -  — ~ 

*  tk>  n+1  ^ 

x,  M  (DZ)K. 

X,M 

0.5  *  [(DZ)k  +  (DZ)K1] 


[Note  that  T^  ^  1  is  a  known  quantity.] 
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Derivation  of  the  Specific  Humidity  Equation 

Again,  neglecting  the  additive  term,  (dq/dt)c>  the  differential  equation  governing 
the  specific  humidity,  q,  may  be  replaced  by  the  difference  equations 


For  k  =  1,  i  =  \;  m  =  fi: 


l,n+l  l,n 

q\,M  ~  q*,M 

(DT) 


fX<q>k,M  '  IY(q)1X,M 


2*0  * 
+  -  * 


(DZ), 


- 

2,  n+1  l,n+l 

- 

K.2  * 

qA.M  qA,iL 

-  VF, 

X,  M 

(°z)2 

X,  [i 

For  k  =  2,  ...,  K  -  1;  £  =  \;  m  =  /i: 
k,  n+1  k,  n 

-  p«<„ 


-  w 


X,  M 


k+l,n+l  k-l,n+l 

rqx.,.  -  V>  1 

(DZ)k+1  +  (DZ)l  J 


2.0 


«DZW  +  <DZ>k] 


k+1,  n+1  k,  n+1 

r  n  —  n 


K 


k+1 
\H  L 


(DZ) 


k+1 


K,  * 
X,H 


k,n+l  k-l,n+l  q 
qA.,M _ V  M 

(DZ),, 


Combining  the  coefficients  of  the  unknowns,  q^,n+1  (k  =  1 . K  -  1),  one  may  put 

A,  n 

the  system  of  difference  equations  in  the  form 


k  k-l,n+l  ,  k  k, n+1  k  k+1, n+1  k 
a.  +  b.  q.  +  c.  q_  =  e.  , 

\,H  X,H  X,n  \,n  X,n  X,n  X,n 

(k  =  1,  ....  K  -  1) 

The  coefficients  a,  b,  and  c  arc  defined  as  for  the  temperature  equation.  The 

k 

coefficients  e,  are  defined  as  follows: 

X,H 


\\  -  +  ‘Y(q)1x!, 


£,  n 
(DT) 


2.0  VF 
(DZ) 


hJL 

2 
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For  k  =  2,  ....  K  -  2: 


\\  =  ‘X(qC  +  IY<q)kkf, 


k,  n 
(DT) 


For  k  =  K  -  1 


K-l 


K-l  .  K, n+1 
w  *  q 

+  +  (DXZ)  ,(DZ)  ' 


2.0  *  K*  *  q^’n+1 

(DZ)  *  [(DZ)  +  (DZ)  ] 


K  '  K-l 
K-l,  n 


(DT) 


Derivation  of  the  Specific  Moisture  Equation 

The  differential  equation  governing  the  specific  moisture  r  is  replaced  by  the 
system  of  difference  equations. 


For  k  =  1,  £  =  X,  m  =  p: 


j.i.n+i  _  rl,n 

=  -  |X(r)l^  •  lY(r)C 


2.0  * 
+  -  * 


(DZ), 

For  k  =  2,  ...,  K  -  1;  £  =  X;  m  =  p: 


k,  n+1  _  k,n 
^  (DT 


1 

2,  n+1  l,n+l 

1 

K?  * 

'rX,p  "  rX,p 

-  VF, 

<DZ>2  J 

X,M 

X,  n 

-  w 


k+l,n+l  k-l, n+1 
fA,p  ~  r\,p 

<DZW  -  <DZ)k 


2.0 


i<DZ>k+i +  (DZ>k' 


K 


k+1 


k+l,n+l  k,n+l 


M  . 


(DZ) 


k+1 


-K,k  ♦ 
\,n 


k, n+1  k-l,  n+1  -| 
"r\.M  ~  r\,/x 


(DZ), 
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Combining  the  coefficients  of  the  unknowns  r^’  n  1  (k  =  1,  K  -  1),  the  system 

X,  p 

of  difference  equations  may  be  put  in  the  form 

k  k-l,n+l  ,  k  k,n+l  k  k+l,n+l  k 
ar  +br  +cr  =e 

X,p  X,p  X,p  X,p  \,n  \,n  h\,n 


(k  =  1,  ....  K  -  1) 

The  coefficients  a,  b,  and  c  are  identical  to  those  defined  above.  The  coefficients  g  are 
defined  as  follows: 


For  k  =  1: 


h!„  -  ix<r"^+  iY(r)C, 


For  k  =  2,  K  -  2: 


rx’n 
X.P  J 

(DT) 
k,  n 

k  k  k  X,  y 

hff,  ‘  +  |Y(r\„  '  pcy 


2.0  VF 


X.  P 


(DZ), 


For  k  =  K  -  1: 


K-l  *  K,  n+1 
w  *  r 

gX,p  1  (  )JX,p  1  (  X,  p  (DZ)  +  (DZ) 


2.0  *  K 


K 


K,  n+1 


hiJk 


(dz)k  *  udz)k  +  (DZ)k  ij 


K  '  K-l 

K-l,  n 

r, 

_ 


(DT) 


Modify  Solutions  for  Phase  Changes 

Application  of  the  method  outlined  above  will  yield  values  for  T,  r,  and  q  at 

t  =  t  .  These  values  are  now  subject  to  an  examination  to  determine  if  they  should 
n+1 

be  modified  due  to  phase  changes  in  the  water  substance.  We  will  repeat  the  technique 
here  to  insure  compatibility  in  our  notation. 


For  k  =  1,  ...,  K;  L  =  1,  ...,  L  and  m  =  1,  ....  M: 


TG  =  T 


,k>n+l  nr  -  k»n+1  nr  -  „k’n+1 

f,  m  ’  QG  -  \m  >  RG  "  r£, m 


OS _ 3,8.*  10;  3 _ 

1.013  -  1.065  •  10-6  *  [Zk  +Ei)m] 


exp 


17.25  * 


TG  -  273. o] 

TG  -  35.7 


CT  = 


[QG  -  QS1  *  ITG  **  21 


5.394  •  103  *  QS  +  4.01  •  10“4  *  [TG  **  2] 


CQ  =  -  4.01  •  10"4  *  CT 
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NON-NEGATIVE 


One  next  performs  the  logic  and  computations  shown  in  Fig.  26. 


Fig.  26.  Logic  and  computations. 


Print  Results 

In  the  section  on  output,  we  indicate  the  character  of  the  output  which  is  desired. 
This  point  in  the  computational  procedure  appears  to  be  a  good  point  at  which  to  print 
the  results  obtained  up  to  t  =  t  . 
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NON-NEGATIVE 


Surface  Temperature  at  t  =  t 

The  value  of  the  air  temperature  at  the  level  of  the  instrument 
shelter  (denoted  by  TS)  will  be  permitted  to  change  in  response  to  diurnal 
radiation  and  to  thermal  advection.  The  computation  is  carried  out  as 
follows. 


For  4=1,  L  and  m  =  1,  ....  M: 

(Not  performed  for  grid  points  over  the  ocean  for  which  =  0.) 


TLAT 

P 

It 


E 

XMR 

TIMER 

If  TIMER  > 

TX 

XX 


If  0  ^  TX  <  2.0  * 


TANF  (SLAT.  ) 

Jt }  m 

SQRTF  [1.  -  (C3  *  (TLAT^  m)  **  2] 

6.0  +  12.0  *  [SQRTF  (1.0  -  P)]  *  [1.5707288  -  0.2121144  * 
P  +  0.0742610  *  (P**2)  -  0.0187293  *  (P**3)]  /  [3.14159] 

TIME/3600.0 

(SLNG  -  75.)/ 15. 

4,  m 

E  -  XMR 

24  set  TIMER  =  TIMER  -  24 
TIMER  -  R 
0. 

(12.0  -  R),  set 


XX  =  1 

ST  =  XX  *  [Cl  sin  (SLAT,  )  -  C2  *  cos  (SLAT.  ) 

4,  m  4,  m 

*  cos  (3.14159  *  TIMER/12.0)] 


YY  =  0 


If  0  s  TX  5  12.0  -  R,  set 
YY  =  1 

RT  =  YY  *  3.14159  *  C2  *  cos  (SLAT.  )  *  sin  (3.14159  *  TIMER/12.0) 

'  4,  m 

II  =  ICLU? 

4,  m 
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For  k  =  1,  ....  K: 

QS  = 


3.8  X  10 


-3 


1.013  -  1.065  X  10“  (Zn  +  Z.  +  E,  ) 

'  n*  —  ■-  £,  m 


exp 


-  273.0, 

17.25  *  |  -y®-: - 

1  T*'n+1  -  35.7 

£,  m 


If  q^’  n  ^  -r  QS  >  0.8  for  any  k,  set 
£,m 


II  =  3 

RADCHC  -  a-Sf  MDT) 
3600 


-1.8  *  RBn  *  TS 

£,  m  I,  m 


+  459.67  *  RB*1  +  RA*1  +  RC,11  *  ST 

£,  m  £,  m  £,  m 


II 


+  RD,  *  RT 
£,  m 


ADVCHG  =  (DT)  *  RATIO.  *  [  -  X(TS),1  -  Y(TS),1  -  w*  *  y] 

'  £,  m  '£,  m  '£,  m  £,  m 


TS.  =  TS,  +  RADCHG  +  ADVCHG 
£,  m  £,  m 


Surface  Relative  Humidity  at  t  =  t 


n+1 


The  value  of  the  relative  humidity  at  the  level  of  the  instrument  shelter  is  computed 
as  follows. 


For  £  =  1,  ....  L;  m  =  1,  ...,  M: 


QS  = 


3.8  X  10 


-3 


[1.013  -  1.065  X  10  (Zn  +  Z.  +  E,  )] 

°£.m  i 


exp 


£,  m 
-  273. 


17.25  * 


r  TS  -  273. 

f Li n _ 1 

l  TS,  -  35.7  j 

It,  m 


1,  n+1 


R®  -  [a  -  CWW)  >  *q4;m  /QS]t(WKT) 


m 
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H  H 

Compute  u  ,  v  ,  and  ICLU  at  t  =  t  , 

_ -  .  -g Kl _ n+1 

H  H 

The  values  of  u^  and  v^  at  t  =  tn  ^  are  to  be  computed  by  quadratic  interpolation 
between  input  quantities  given  at  six-hour  intervals.  Let  us  denote  these  input  quantities 
by  UGHljj)  m;  UGH2£  m;  UGH3£m;  VGHl£m;  VGH2£  m;  VGH3£m.  The  computation  may 
then  proceed  as  follows:  set 

TAU  =  (n  +  1)  *  (DT)/(6.0  *  3600.0) 

For  £  =  1,  ...,  L;  and  m  =  1 . M: 

AF  =  [4.  *  UGH2  -  UGH3.  -  3.  *  UGH1„  ]  /  2.0 
£,  m  £,  m  £,m 

BF  =  [UGH3.  -  2.0  *  UGH2.  +  UGII1,  ]  /2.0 

£,  m  £,m  £,  m 

H 

U  =  UGH1  +  AF  *  TAU  +  BF  *  (TAU  **  2) 

g  £,  m  £,  m  '  ' 

AF  =  [4.  *  VGH2.  -  VGH3  -  3.  *  VGH1  ]  ,  .  A 
£,  m  £,  m  £,m  /  2.0 

BF  =  [VGH3.  -  2.0  *  VGH2  +  VGH1.  ]  /  2.0 
£,  m  £,m  £,  m 

vH,  =  VGH1„  +  AF  *  TAU  +  BF  *  (TAU  **  2) 
g£,m  £,m 

To  compute  the  appropriate  value  of  ICLU  at  t  =  t  we  make  use  of  character¬ 
istic  values  of  ICLU  over  six-hour  intervals.  Let  these  input  quantities  be  denoted  by 


ICLU1.  and  ICLU2„  . 

£,  m  £,  m 

Having  computed  TAU  as  above,  ICLU^  is  determined  as  follows: 
If  TAU  si,  set 

ICLU.  =  ICLU1 


£,  m 

If  TAU  >1,  set 


£,  m 


ICLU .  =  ICLU2 

£,  m  £,  m 


Replace  Values  at  t^  by  Values  at  t^+^ 

The  newly-computed  time-dependent  quantities  may  now  replace  the  values  pre¬ 
viously  stored  at  level  n.  The  old  values  are  lost.  The  locations  reserved  for  the 
unknowns  at  n  +  1  should  now  be  set  equal  to  zero. 
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Step  Time 


Set  TIME  =  TIME  +  DT  and  adjust  indices  as  required. 
(Check  if  Computation  is  Completed) 
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